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1 . 0  INTRODUCTION 


The  Electrothermal  Gun  (ETG)  offers  a  way  to  circumvent  a  fundamental 
performance  limitation  associated  with  conventional  chemical  propulsion. 
Since  the  specific  chemical  energy  of  conventional  propellants  is  limited, 
increases  in  total  propulsive  energy  are  accompanied  by  a  concomitant 
increase  in  the  total  mass  of  the  propellant.  As  the  propellant  is  required 
to  follow  the  projectile,  a  fraction  of  the  total  energy  is  necessarily 
scored  as  kinetic  energy  of  the  moving  propellant.  This  kinetic  energy 
represents  a  loss  of  efficiency  which  is  proportional  to  C/M,  the  ratio  of 
the  mass  of  the  propellant  to  that  of  the  projectile.  For  conventional 
ammunition  in  which  the  muzzle  velocity  is  approximately  1  km/s ,  the  value 
of  C/M  is  about  0.25  and  the  kinetic  energy  loss  is  not  large.  However,  for 
hypervelocity  weapons  operating  at  muzzle  velocities  of  the  order  of  3  km/s, 
the  value  of  C/M  can  be  significantly  greater  than  1  at  which  point  the 
kinetic  energy  losses  can  begin  to  dominate:  a  regime  of  diminishing  returns 
is  entered  as  the  chemical  energy  is  expended  primarily  to  accelerate  the 
propellant  itself  rather  chan  the  projectile. 

The  ETG  uses  an  intensely  heated  plasma  which  is  injected  into  a  work¬ 
ing  fluid  to  create  a  pressurized  gas  which  in  turn  propels  the  projectile. 
The  basic  configuration  is  illustrated  in  Figure  1.1.  The  review  by  Juhasz 
et  al  [1]  may  be  consulted  for  a  discussion  of  the  various  types  of  config¬ 
urations  presently  under  consideration;  these  include  examples  in  which  the 
working  fluid  is  either  inert  or  reactive.  As  shown  in  Figure  1.2,  the 
working  fluid  may  have  a  central  cavity  and  Che  reactive  examples  may  util¬ 
ize  richer  monopropellants  or  bipropellants.  Since  the  plasma  is  created  by 
an  external  supply  of  electrical  energy,  it  is  possible,  in  principle,  to 
achieve  an  arbitrarily  high  energy  density  in  the  heated  working  fluid. 
Accordingly,  the  C/M  constraint  associated  with  conventional  chemical  pro¬ 
pulsion  is  removed.  However,  a  new  constraint  appears  in  connection  with 
the  thermal  response  of  the  Cube.  It  is  necessary  Chat  the  heat  capacity  of 
the  working  fluid  be  sufficiently  high  as  to  moderate  the  temperature  of 
working  fluid  so  that  the  heat  transfer  to  the  gun  tube  is  not  excessive. 
Since  Che  heat  capacity  is  obviously  proportional  to  the  total  mass  of  the 
heated  fluid,  it  is  clear  that  this  thermal  constraint  is  in  effect  a  C/M 
constraint  similar  to  that  associated  with  the  kinetic  energy  loss.  How¬ 
ever,  the  ETG  admits  the  use  of  working  fluids  having  much  larger  specific 
heat  capacities  chan  those  of  chemical  propellants.  It  is  this  additional 
degree  of  freedom  of  design  which  preserves  the  apparent  propulsive  advan¬ 
tage  attached  to  the  use  of  an  external  supply  of  energy.  Otherwise,  keep¬ 
ing  the  heat  transfer  to  the  tube  comparable  to  that  obtained  with  conven¬ 
tional  propellants  would  require  a  mass  of  working  fluid  which  diluted  the 
plasma  energy  to  a  value  comparable  to  that  of  the  conventional  propellant 
itself,  at  which  point  the  conventional  level  of  kinetic  energy  loss  would 
be  experienced. 


1.  Juhasz.  A.,  Jamison,  K. ,  White,  Y.  and  Wren,  G.  "Introduction  to 
Electrothermal  Gun  Propulsion"  Proceedings  of  the  25th  Jannaf 
Combustion  Meeting  1988 
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(a)  Single  Flxiid.  Homogeneous 


(b)  Single  Fluid,  with  Cavity 


(c)  Bipropellant,  no  Cavity 


(d)  Bipropellant,  with  Cavity 


Figure  1.2  Types  of  ETC  Configurations 
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The  preceding  discussion  has  been  based  on  an  implicit  assumption, 
namely  that  the  plasma  energy  is  uniformly  distributed  throughout  the 
working  fluid.  In  fact,  this  is  not  likely  to  be  true.  Thus  one  must  ask 
to  what  extent  are  the  ballistic  performance  and  the  rate  of  heat  transfer 
to  the  tube  affected  by  the  details  of  the  mixing  process.  It  may  be 
anticipated  that  propulsion  will  be  affected  less  strongly  by  the  details  of 
mixing  than  will  the  heat  transfer.  Since  the  spacemean  pressure  is  essent¬ 
ially  proportional  to  the  total  energy  within  the  mixing  chamber,  the  over¬ 
all  rate  of  pressurization  is  not  expected  to  depend  strongly  on  the  rate  of 
mixing.  Such  dependence  as  arises  will  stem  from  non- idealities  in  the 
equation  of  state,  the  influence  of  the  effective  C/M  on  the  pressure  grad¬ 
ient  in  the  gun  tube,  and  on  transient  wave  phenomena.  However,  the  mixing 
process  will  be  of  paramount  importance  in  regard  to  the  distribution  of 
temperature  and  therefore  the  heat  transfer  to  the  tube.  Poor  mixing  will 
produce  intensely  heated  gas  but  at  the  same  time  will  preserve  a  region  of 
relatively  cool  fluid.  Heat  transfer  to  the  tube  may  be  either  reduced  or 
increased  by  poor  mixing  depending  on  whether  or  not  the  cool  fluid  acts  as 
an  insulating  layer  near  the  tube  wall. 

In  Chapters  2.0  and  3.0  of  this  report  we  focus  attention  on  the 
ballistic  consequences  of  the  rate  of  mixing  insofar  as  this  topic  can  be 
addressed  in  terms  of  lumped  parameter  and  one -dimensional  continuum  models. 
Meaningful  analysis  of  the  implications  of  the  mixing  rate  on  the  heat 
transfer  to  the  tube  requires  a  two-dimensional  model.  Chapter  4.0  presents 
the  formulation  of  such  a  model. 

The  results  presented  in  Chapters  2.0  and  3.0  are  all  based  on  the 
assumption  that  the  rate  of  injection  of  the  plasma,  and  its  energy  content, 
are  predetermined.  Only  the  state  of  the  mixture  of  the  plasma  and  the 
working  fluid  is  modeled  explicitly.  The  approach  is  therefore  similar  in 
spirit  to  that  described  by  Oberle  [2]  who  used  the  lumped  parameter  inter¬ 
ior  ballistics  code  IBHVG2  [3]  to  model  the  response  of  the  working  fluid  to 


2.  Oberle,  W.  "Electrothermal  Guns  --  A  Theoretical  Investigation 

of  Factors  for  Optimal  Performance"  Ballistic  Research  Laboratory 
Report  BRL-TR-2999  1989 

3.  Anderson,  R.D.  and  Fickle,  K.D.  "IBHVG2  -  A  User's  Guide"  BRL 

Technical  Report  BRL-TR-2829  Ballistic  Research  Laboratory,  Aberdeen 
Proving  Ground,  HD  1987 
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Che  addition  of  the  plasma.  Models  which  address  specifically  the  formation 
of  the  plasma  have  been  advanced  by  Powell  and  Zielinski  [4]  ,  Tidman  eC  al 
[5],  Loeb  and  Kaplan  [6]  and  by  Chryssomallis  et  al  [7]. 

Chapter  2.0  contains  a  complete  discussion  of  the  lumped  parameter 
model  developed  under  the  present  contract  Cask.  The  model  is  encoded  as  a 
program  referred  to  as  LUMPET.  Appendix  A  contains  a  listing  of  the  code 
together  with  a  complete  description  of  the  input  files.  Sample  data  bases 
and  solutions  are  presented  in  Appendices  C  and  D.  In  addition  to  the 
discussion  of  the  model.  Chapter  2.0  presents  a  number  of  numerical  results. 
The  model  is  used  first  to  predict  the  influence  of  the  rate  of  mixing  on 
the  interior  ballistics  for  a  nominal  ETC  configuration.  Second,  the  code 
is  applied  in  an  inverse  mode  to  deduce  the  rate  of  mixing,  or  decomposition 
of  the  working  fluid,  when  both  the  plasma  flux  and  the  gun  pressure  history 
are  given.  The  first  set  of  results  shows  relatively  little  sensitivfty  of 
the  interior  ballistics  Co  the  rate  of  mixing.  Not  surprisingly,  the  second 
set  of  results  make  it  clear  that  very  accurate  experimental  data  will  be 
required  if  the  code  is  to  be  used  in  an  inverse  mode  to  determine  the  rate 
of  mixing  with  any  degree  of  precision. 

To  clarify  further  the  dependence  of  pressure  on  the  extent  of  decom¬ 
position  of  the  working  fluid  we  provide  a  short  computer  code  referred  to 
as  PMAP.  The  code  is  described  in  Section  2.3  and  a  listing,  together  with 
a  description  of  input,  is  given  in  Appendix  B. 

The  discussion  of  Chapter  3.0  provides  only  a  brief  summary  of  the  one- 
dimensional  model.  Reference  is  made  Co  independent  documentation  of  the 
relevant  code,  referred  to  as  XNOVAKTC  or,  more  briefly,  as  XKTC  [8]. 
Chapter  3.0  only  presents  those  model  details  which  are  particular  to  the 
representation  of  the  ETC  and  which  were  encoded  as  part  of  the  present 
contract  task.  The  revised  version  of  XKTC  is  used  to  assess  the  influence 
of  Che  rate  of  mixing  from  a  one -dimensional  perspective  which  allows 
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spatial  non- uniformity.  We  study  the  same  ETG  configuration  as  in  Chapter 
2.0.  The  results  confirm  those  of  Chapter  2.0  in  the  sense  that  only  modest 
sensitivity  is  seen  for  the  sample  problem  under  consideration.  A  sample 
XKTC  data  base  is  presented  in  Appendix  E. 

In  Chapter  4.0  we  present  the  details  of  a  two-dimensional  model  of  the 
ETG.  The  discussion  includes  a  statement  of  the  principal  governing 
equations  and  of  the  intended  isethod  of  solution.  However,  the  relevant 
coding  and  the  determination  of  ntimerical  solutions  is  outside  the  scope  of 
the  present  task  and  will  be  presented  in  a  subsequent  report. 
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2.0  LUMPED  PARAMETER  MODEL 

This  chapter  has  three  sections.  In  Section  2.1  we  present  the  details 
of  a  lumped  parameter  model  of  the  interior  ballistics  of  the  ETC.  In 
Section  2.2  we  present  some  numerical  results  which  illustrate  the  depend¬ 
ence  of  the  interior  ballistics  on  the  rate  of  mixing  of  the  plasma  with  the 
working  fluid.  In  Section  2.3  we  consider  the  use  of  the  model  as  a  means 
of  determining  the  rate  of  mixing  when  the  pressure  history  of  the  gun  is 
given  in  addition  to  a  characterization  of  the  plasma  flux.  As  we  shall  see 
in  Section  2.2,  the  interior  ballistics  of  the  ETC  are  only  weakly  dependent 
on  the  rate  of  mixing.  This  is  a  desirable  result  from  the  point  of  view  of 
the  gun  designer  since  it  implies  that  pressurization  is  dominated  by  the 
rate  of  addition  of  plasma.  On  the  other  hand,  it  implies  that  the  deter¬ 
mination  of  the  rate  of  mixing  from  the  pressure  history  will  require  rela¬ 
tively  precise  data.  We  explore  thia  latter  point  in  some  detail  in  Section 
2.3. 


2 . 1  GOVERNING  EQUATIONS 

We  assume  that  we  have  a  mixing  chamber  where  volume,  prior  to  motion 
of  the  projectile,  is  V^.  The  chamber  contains  a  total  mass  C  of  working 
fluid  and  an  initial  mass  m^  of  gas.  The  density  of  the  unmixed  portion  of 
the  working  fluid  is  Pp(p),  a  function  of  the  pressure,  p.  At  any  time  t,  a 
fraction  a(t)  of  the  working  fluid  has  been  vaporized  and  mixed  uniformly 
with  the  injected  plasma.  While  we  neglect  the  time  delay  associated  with 
the  mechanical  process  of  mixing  each  element  of  vaporized  working  fluid 
with  the  mixture  of  gases  already  present,  we  do  consider  the  possibility 
that  the  approach  to  chemical  equilibrium  requires  a  period  of  time  which  is 
not  negligible  from  an  interior  ballistic  perspective.  Accordingly,  the 
mixture  of  gases  is  formally  decomposed  into  an  intermediate  species  Sj, 
whose  mass  fraction  is  Y,  and  a  final  species  Sp,  whose  mass  fraction  is 
obviously  1  -  Y.  At  any  given  time,  the  total  mass  of  plasma  in  the  mixing 
chamber  is  m^(t)  and  the  energy  of  the  plasma  is  E^Ct). 

As  noted  previously,  we  assume  that  m^Ct)  and  Ej^(t)  are  furnished  as 
tabular  data;  the  plasma  flux  is  not  modeled  directly.  In  subsequent 
subsections  we  discuss  the  equations  of  motion  of  the  projectile,  the 
relations  governing  the  pressure  gradient,  the  equations  of  state  of  the 
unvaporized  working  fluid  and,  finally,  the  equations  which  govern  the 
mixture  of  gases. 

Motion  of  the  Projectile 

The  projectile  has  mass  M,  displacement  x^  and  velocity  v  .  Its  motion 
is  resisted  by  two  forces,  the  first  is  due  ^to  the  engraving  band  and  is 
represented  by  a  pressure  Pj._g  while  the  second  is  due  to  the  compression  of 
air  in  front  of  the  projectile  and  is  represented  by  The  equation  of 
motion  of  the  projectile  is 


dv  g  A. 

p  _  ^o  D  .  _ 

dt  M  ^^base  ^res  ^air^ 


(2.1) 
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where  is  the  bore  area  and  Is  Che  pressure  at  the  base  of  Che 
projectile.  We  assume  that  Pfag  ~  Pres^^^  Chat  Che  functional  depend¬ 
ence  Is  characterized  by  a  table  of  values.  The  value  of  Pg^j.  Is  assumed 
Co  follow  from  the  jiimp  conditions  for  a  shock  wave  In  which  the  gas  ahead 
of  Che  shock  Is  at  rest  trhlle  the  gas  behind  the  shock  has  a  velocity  equal 
Co  that  of  the  projectile.  It  follows  that  [9] 


Pair  -  Pa 


(1  +  Mg‘) 


p  p  a  a 


1  ‘/* 


2(1 


2. 

U  )c 
a  I 


-  M. 


(2.2) 


where  -  (y^  -  l)/(y  +  1)  and  i-s  the  ratio  of  specific  heats  of  the 
air  In  the  barrel.  We  ^so  have  c^,  the  speed  of  sound  In  the  unshocked  air 
and  p  Is  Che  pressure.  The  tiser  of  the  code  has  the  option  of  Including 
the  effect  of  Che  air  shock  or  of  neglecting  It.  The  properties  of  Che  air 
are  assumed  to  be  characterized  by  the  values  y.  •  1.4,  P-  *■  0.101  MPa  and 
c^  -  343.2  m/s. 

Pressure  Gradient  Relations 

The  eqixatlon  of  motion  of  Che  projectile  Involves  Che  base  pressure 
Pbgge*  balance  equations  which  describe  the  thermodynamic  state  of  the 
mixture  will  be  predicated  on  Che  spacemean  pressure,  p.  Moreover,  It  Is  of 
Interest  Co  the  charge  designer  to  know  the  breech  pressure,  p^_,  since  this 
Is  expected  to  be  Che  highest  pressure  In  the  gun.  The  model  admits  two 
different  sets  of  analytical  relations  between  these  values  of  pressure. 
The  first  set  Is  referred  Co  as  the  Lagrange  relations  [10]  while  the  second 
set  Is  referred  to  as  the  Chambrage  relations  [11].  The  first  set  Is  very 
familiar  to  Interior  balllstlclans  and  cakes  no  account  of  the  variations  In 
cross-sectional  area  of  the  chamber.  The  second  set  reflects  the  non- 
uniform  area  of  the  chamber.  Relations  of  Che  second  type  appear  to  have 


9.  CouranC,  R.  and  Friedrichs,  K.O.  "Supersonic  Flow  and  Shock  Waves" 

Interscience ,  New  York  1948 

10.  Comer,  J.  "Theory  of  the  Interior  Ballistics  of  Guns" 

Wiley,  New  York  1950 
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been  reported  first  by  Vintl  [12].  Slnllar  results  have  been  obtained  nore 
recently  by  Morrison  and  Wren  [13]  and  cast  Into  a  more  general  format  by 
Robbins  et  al  [11]. 

Both  the  Lagrange  and  the  Chambrage  relations  are  based  on  the  assump¬ 
tion  that  the  flow  Is  homogeneous.  However,  In  the  present  application  we 
cannot  assume  that  the  unvaporlzed  working  fluid  Is  necessarily  at  rest. 
Nor  can  we  assume  that  it  moves  with  the  same  velocity  as  the  mixture  of 
gases.  We  therefore  Introduce  a  parameter  fi(t) ,  referred  to  as  the  Lagrange 
coupling  coefficient,  which  reflects  the  degree  to  which  the  unvaporlzed 
working  fluid  moves  with  the  gases.  If  ^  -  0,  the  vinvaporlzed  working  fluid 
Is  assumed  to  be  at  rest.  If  ^  -  1,  the  working  fluid  moves  with  the  gases. 
The  value  of  /S  Is  required  to  be  specified  by  the  user  of  the  code.  In 
Section  2.3  we  discuss  the  experimental  determination  of  p(t)  from  an 
Inverse  analysis  of  the  observed  Interior  ballistic  behavior. 

Lagrange  Relations 

Let  K  be  defined  as 


K 


C(a  +^(1  -  a))  + 
M 


(2.3) 


The  value  of  the  breech  pressure  and  the  base  pressure  are  related  to  the 
spacemean  pressure  according  to 


_  1 

^base  "  1  +  K/3  ^ 


(2.4) 


Pbr 


1  +  K/2 
1  +  K/3  P 


(2.5) 


In  addition,  we  have  the  following  expression  for  the  kinetic  energy  of  the 
unvaporlzed  working  fluid  and  the  mixture  of  gases, 


^in 


KM  2 

6go  P 


(2.6) 
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Chanbrage  Relations 

Again,  with  K  as  defined  by  Equation  2.3  we  have  the  following  rela¬ 
tions  bet  .en  the  various  values  of  pressure. 


^base  k 

2 


(2.7) 


Pbr 


k  k 

3  X 

k 


2 


+  k 

4 


(2.8) 


The  kinetic  energy  is  given  by 


^in 


v’(Zp) 


(2.9) 


where  z.  is  the  distance  from  the  breech  face  to  the  base  of  the  projectile 
and  we  Introduce  the  following  expressions. 


V(z) 


(2.10) 


Jj(z) 


fZ 

dz 
o 


(2.11) 


J^(z)  -  V*(z)/A*(z) 


(2.12) 


J3(z) 


I  A(z)Jj(z)dz 
*  o 


(2.13) 
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f^(z)  -  I  A(z)Jj(z)d2 


KMA. 

a  (t) - ^ 

»  ..2 


8o^  <»p>  >-  P 


Vo  8o^» 


V(z^)  *  M  ^Pres  *  ^air^ 


VO  -  - 


V‘(Zp) 


b(t)  -  - 


A  2  2 

KM  \  \ 
28o 


a  J,  bJ 

k.  -  a  J.  +  bJ  *  ^  - * 


I  »  I  --2  V(Zp)  V(Zp) 


k  -  1  - 


a  J 

2  3 

a  J  +  — — 7 
2  I  V(z  ) 


k  -  1  - 


a  J 
2  1 


k  -  -  a  J  -  bJ 

4  II  2 


(2.14) 


(2.15) 


(2.16) 


(2.17) 


(2.18) 


(2.19) 


(2.20) 


(2.21) 
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Assuaing  Chat  Che  bulk  modulus,  Kp,  of  Che  unvaporlzed  working  fluid  is 
at  most  a  linear  function  of  pressure 


K  -  K  +  K  (p  -  p  ) 
p  t  2'*^  *^o 


(2.22) 


we  have  Che  equation  of  state  in  the  fora. 


P(P)  - 


(1  +  (P  -  Po)/K,) 
Pft  ^  ^ 


^  a*/  (V  - 

*^o  i 


if  K  -  0 
2 


(2.23) 


if  K  #0 
2 


where  p  -  Pn(Po)  •  Corresponding  Co  this  equation  of  state  we  have  Che 
Pq  P  o 

following  expression  for  the  energy  per  unit  mass  stored  by  the  compressed 
fluid, 


ep(p)  - 


<P  ’  Po> 

Pp(P) 


if  Kg  -  0  or  1, 


(2.24) 


(p  -  p  )  K  +  K  (p  -  p  )  K 

Pp(P)  (Kg  -  l)Pp(p)  ■  (K^  -  l)Pp^ 


if  Kg  #  0 
and  Kg  #  1. 


Equations  for  the  Mixture  of  Plasma  and  Unvaporized  Working  Fluid 

At  any  given  time  the  mixture  of  plasma  and  vaporized  working  fluid 
occupy  Che  instantaneous  volume 


V 


V 

o 


+  X 

P 


\  -  (1  -*>  f" 

^  Pp 


(2.25) 
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The  density  of  the  mixture  Is 


P 


m.  +  a  +  atC 
1  o _ 

V 


(2.26) 


where  p  Is  understood  to  be  a  spaceaean  value  like  the  pressure  p. 

We  may  write  a  balance  equation  for  the  Intermediate  species  Sj  In  the 

form 


dt 


(pYV) 


-  VAT"(Yp)''exp 


(2.27) 


Evidently,  the  first  term  on  the  right  hand  side  of  Equation  2.27  repre¬ 
sents  production  of  Sj  via  decomposition  of  the  working  fluid  while  the 
second  term  represents  depletion  according  to  an  Arrhenius  reaction  law. 
The  reaction  term  Introduces  A,  the  pre -exponential  factor;  n,  the  temper¬ 
ature  exponent;  v,  the  reaction  order;  E»,  the  activation  energy;  and  R^, 
the  universal  gas  constant.  Using  the  balance  of  mass  for  the  mixture  as  a 
whole,  Equation  2.27  can  be  recast  Into  the  computational  form 


dt 


(1  -  Y)  ^  Y  ^1 
pV  ^  dt  ■  pV  dt 


1  ._n,„  .V 
-  -  AT  (Yp)  exp 


(2.28) 


The  specific  heats,  covolume  and  molecular  weight  of  the  mixture  can  be 
written  as  follows 


c 

V 


+  (1  -  Y)c 

V 


F 


c 

P 


+  (1  -  Y)c^ 

pf 


b  -  Ybj  +  (1  -  Y)bp 


M 


w 


1 

+  (1  -  Y) 

M 


(2.29) 


(2.30) 


(2.31) 


(2.32) 
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Here  and  c  are  respectively  the  specific  heats  at  constant  volume  and 
constant  pressure,  b  is  the  covolume  and  is  the  molecular  weight.  Ve 
have  used  subscripts  I  and  F  to  denote  properties  of  and  Sp  respectively. 
It  should  be  noted  that  we  neglect  the  influence  of  m^,  the  mass  of  gas 
initially  present.  The  properties  of  Sp  are  understood  to  reflect  the 
mixture  of  plasma  and  working  fluid  at  chemical  equilibrium. 

We  assume  that  the  spacemean  properties  satisfy  the  covolume  equation 
of  state  in  the  form 


P  - 


(2.33) 


where  is  the  effective  internal  energy  whose  determination  we  now 

discuss . 

We  may  write  a  global  energy  balance  for  the  mixture  as  follows 


[m.  +  +  aCJId  -  YXe^  +  Hp)  +  Y(ej  +  Hj)]  +  (1  -  a)Cep(p) 

"  2i;  V  *  "kin  *  ^  •  »>“o 

-  +  Cep(p^)  +  CH^  +  +  “o)t(l  -  Y)Hp  +  YH^] 


Here  Cp  and  Hp  are  the  thermal  energy  and  the  heat  of  formation  of  the  final 
products  Sp,  and  ex  and  Hp  are  the  corresponding  properties  of  the  intermed¬ 
iate  products.  The  second  term  on  the  left  hand  side  is  the  energy  stored 
in  the  unvaporlzed  working  fluid.  The  third  term  is  the  kinetic  energy  of 
the  projectile  and  the  fourth  is  the  kinetic  energy  of  the  mixture  of  gases 
and  the  unvaporlzed  working  fluid.  The  fifth  term  is  the  work  done  against 
the  forces  resisting  the  motion  of  the  projectile. 


W 


f 


X 

P 


‘b'^res  *  Pair*'*' 


o 


(2.35) 
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The  sixth  term  represents  the  heat  lost  to  the  tiibe  and  is  presently  set 
equal  to  zero.  The  last  term  on  the  left  hand  side  represents  the  heat  of 
formation  of  the  unvaporized  working  fluid.  The  terms  on  the  right  hand 
side  correspond  respectively  to  the  initial  energy  of  the  gas  in  the  coo^us- 
tion  chamber,  the  energy  added  by  the  plasma,  the  energy  of  compression  of 
the  working  fluid  at  the  initial  conditions,  the  heat  of  formation  of  the 
working  fluid  and  finally,  a  heat  of  formation  term  consistent  with  our 
neglect  of  the  detailed  influence  of  the  initial  mass  of  gas  and  the  mass  of 
the  plasma. 

Ue  introduce 


Q„  -  Hj  .  Hj.  ,  (2.36) 


Q,  -  ,  (2.37) 


and  we  write 


e  -  C  T  ,  (2.38) 

I  Vj 

where  T  is  the  temperature  of  the  mixture.  Then  we  solve  for  ep  in  the  form 


ep  -  Ei  +  0aC(Qjj  -  Q^)  -  aCYQj,  - 


-  (m,  +  m  +  oC)YC  T  -  (1  -  o)Ce 
i  o  P 

(2.39) 

In  Equation  2.39  we  have  introduced  a  coefficient  0  whose  value  depends  on 
the  convention  adopted  in  the  determination  of  the  thermodynamic  properties. 

The  model  admits  two  distinct  conventions.  The  first  is  referred  to  as 
the  constant  thermochemistry  option  and  takes  the  properties  of  Sp  to  be 
fixed  at  all  times.  In  this  case,  0-1,  and  the  effective  ground  state  for 
the  internal  energy  is  that  defined  by  the  heat  of  formation  of  the  unvapor¬ 
ized  working  fluid.  Thus  Q.^  is  to  be  interpreted  as  the  heat  required  to 
vaporize  a  unit  mass  of  working  fluid  and  Qu  is  to  be  interpreted  as  the 
heat  released  per  unit  of  mass  transformed  from  to  Sp.  The  effective 
internal  energy  as  required  by  Equation  2.33  is  taken  to  be  equal  to  ep 

as  defined  in  Equation  2.39.  In  the  second  convention,  referred  to  as  the 
variable  thermochemistry  option,  we  set  the  value  of  0  -  0.  In  this  case 


/[(m^  +  m^+  aC)(l  -  Y) ] 


*  \in  -  «f  *  ^  ^  ^o  %(Po> 
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the  effective  ground  state  is  that  of  the  final  products  following  heating 
by  the  plasma.  The  properties  of  the  mixture  are  assumed  to  have  been 
characterized  by  a  series  of  runs  using  the  BIAKE  Code  as  described  by 
Oberle  [2].  A  table  of  values  is  presumed  to  be  given  in  which  the  effec¬ 
tive  thermal  energy,  e^££,  and  the  values  of  y  -  c^/c^,  b,  and  Qy  are 
specified  for  a  series  of  values  of  energy  added  per  unit  mass.  These 
correspond  to  the  thermodynamic  equilibritim  of  the  mixture  and  accordingly 
discount  the  energy  associated  with  vaporization  of  the  working  fluid  and 
reaction  of  any  intermediate  products.  In  this  case  is  interpreted  to  be 
the  difference  in  heat  of  formation  of  the  intermediate  and  final  products 
and  is  positive  if  the  transformation  of  Sj  into  Sp  is  exothermic.  The 
value  of  ep  determined  by  Equation  2.39  is  used  to  determine  eg££  and  the 
other  variables  by  linear  interpolation  of  the  tabulated  BIAKE  Code  results . 

Solution  of  the  Governing  Equations 

The  model  consists  of  a  number  of  algebraic  equations  and  a  total  of 
four  differential  equations,  two  for  the  equations  of  motion  of  the 
projectile,  one  for  the  work  done  against  the  forces  resisting  projectile 
motion  and  one  for  the  rate  of  change  of  the  mass  fraction  of  the 
intermediate  species.  The  differential  equations  are  solved  using  a 
standard  fourth  order  Runge-Kutta  algorithm  [14]  with  a  fixed  time  step. 

Several  of  the  algebraic  relations  are  non-linear.  The  solution  of  the 
energy  equation  is  accomplished  by  means  of  Newton's  method  with  the 
pressure  as  the  independent  variable.  If  a  chemical  reaction  is  considered, 
a  solution  by  Newton's  method  is  also  required  in  order  to  determine  the 
temperature.  The  model  is  encoded  into  a  Fortran  program  referred  to  as 
LIMPET.  Appendix  A  contains  a  listing  of  the  code  together  with  a  discus¬ 
sion  of  the  input  data.  Sample  solutions  are  presented  in  Appendices  C  and 
D.  These  pertain  to  the  discussion  of  Sections  2.2  and  2.3  respectively. 

In  its  normal  mode  of  operation  the  code  accepts  a  tabular  description 
of  the  time  dependence  of  a(t),  m£(t)  and  E^Ct).  From  these  inputs, 
together  with  a  characterization  of  the  thermodynamic  properties  of  the 
mixture,  the  interior  ballistic  history  is  computed.  In  a  second  mode  the 
user  may  specify  a  constant  value  of  breech  pressure  together  with  a  total 
value  of  plasma  mass  and  energy.  The  code  will  then  determine  the  histories 
a,  m£  and  E^  to  provide  the  desired  constant  pressure  up  to  the  point  at 
which  the  plasma  supply  is  exhausted.  The  rate  of  decomposition  is  taken  to 
be  strictly  proportional  to  the  plasma  energy  flux.  A  midpoint  search  is 
used  to  determine  the  value  of  E^  at  each  time  step. 

In  a  third  mode  of  operation  the  code  accepts  a  tabular  description  of 
the  observed  pressure  history  in  the  gun  together  with  the  details  of  m£(t) 
and  E^Ct).  This  is  referred  to  as  the  inverse  mode  and  the  code  proceeds  to 
deduce  ot(t)  and,  if  both  breech  and  base  pressure  are  given,  /S(t)  as  well. 
We  provide  further  details  of  this  mode  of  operation  and  the  method  of 
solution  in  Section  2.3. 


14.  Ralston,  A.  "A  First  Course  in  Numerical  Analysis"  McGraw-Hill  1965 
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2.2  Ballistic  Implications  of  Rate  of  Mixing 


In  Che  present  section  we  illustrate  the  use  of  the  lumped  parameter 
model  Co  probe  the  sensitivity  of  the  interior  ballistics  of  the  ETG  to  the 
rate  of  decomposition  of  the  working  fluid.  We  recall  that  the  model 
assumes  the  working  fluid  to  be  mixed  with  the  ambient  as  soon  as  it 
decomposes.  Thus  the  discussion  of  this  section,  and  the  next,  may  be 
viewed  as  addressing  either  the  rate  of  decomposition  or  the  rate  of  mixing 
since  there  is  no  difference  between  the  two.  The  study  is  similar  to  Chat 
reported  in  a  previous  paper  [15]  in  which  a  rudimentary  version  of  the 
present  model  was  employed.  The  numerical  results  presented  here  differ 
from  chose  of  Reference  IS  due  Co  revisions  to  the  model  and  certain  minor 
differences  in  the  data  bases. 

We  first  discuss  the  influence  of  the  rate  of  sensitivity  in  general 
analytical  terms  by  focussing  on  Che  implications  of  Che  equation  of  state. 
We  then  discuss  the  data  bases  and  the  numerical  results. 

Dependence  of  Pressure  on  a 

The  influence  of  Che  rate  of  mixing  can  be  anticipated  from  an  examina¬ 
tion  of  Equation  2.33  which  we  repeat  as 


ir 


p  - 


i  -  b 

P 


As  discussed  in  the  previous  section,  we  admit  two  conventions  with  respect 
to  Che  representation  of  the  thermochemical  properties  of  Che  mixture  of 
plasma  and  vaporized  working  fluid,  namely  a  constant  thermochemistry  option 
and  a  variable  thermochemistxy  option.  In  the  latter,  the  effective  energy, 
&^ff,  and  Che  values  of  y  and  b  are  functions  of  ep,  which  is  seen  to  be  a 
strong  function  of  a  through  the  denominator  of  Equation  2.39.  In  the 
former  case,  however,  we  have  y  and  b  constant  and  “  ep.  We  therefore 
confine  our  discussion  for  the  moment  to  the  former  case  and  make  the 
simplifying  assumption  chat  Q^,  Q,^,  W^,  Q^,  E^,  Op,  m^,  m^  and  Y  are  all 
zero.  That  is  to  say,  we  assume  the  working  fluid  to  be  inert  and  the 
losses  associated  with  work  against  friction,  heat  loss  Co  Che  tube  and 
compression  of  the  unvaporlzed  working  fluid  to  be  negligible.  We  also 
neglect  the  mass  of  the  plasma  and  the  influence  of  the  initial  ambient. 
With  these  assumptions,  Equations  2.26  and  2.39  show  that 


e 


eff 


aC 


\in 


(2.40) 


15.  Gough,  P.S.  "Influence  on  Interior  Ballistics  of  Electrothermal 
Gun  of  Rate  of  Mixing  of  Plasma  with  Working  Fluid"  Proceedings 
of  the  26th  Jannaf  Combustion  Meeting  1989 
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We  use  Che  siiiple  Lagrange  relation  for  Che  kinetic  energy  as  given  in 

Equation  2.6  whereupon  2.33  and  2.40  give  the  result 


(2.41) 


In  this  simplified  sitxiation  we  have  the  following  conclusions.  First,  p  is 
completely  insensitive  to  a  in  the  special  case  when  —  1  and  1/Pn  ~ 
Second,  p  will  be  only  weakly  sensitive  to  a  in  the  more  general  situation 
defined  by  the  pair  of  inequalities. 


C/H  «  1 


V 

o 


(2.42) 


From  an  examination  of  the  second  inequality  it  is  apparent  that  the 
presence  of  initial  ullage,  >  C/p  ,  may  be  expected  to  desensitize  the 
ETC  interior  ballistics  to  Che  rate  of  mixing.  In  general,  configurations 
having  a  low  C/M  and  low  initial  loading  density  of  the  working  fluid  will 
be  expected  to  exhibit  less  sensitivity  to  the  rate  of  decomposition  than 
configurations  with  large  values  of  C/N  and  high  loading  densities.  Of 
course,  all  this  discussion  is  predicated  on  the  assumption  that  the 
thermochemical  properties  are  constant.  As  we  shall  see  in  the  numerical 
results  which  follow,  the  variation  in  properties  with  variations  in  the 
energy  density  ep  may  exert  a  dominant  role  as  regards  the  sensitivity  of 
the  interior  ballistics  to  the  rate  of  mixing. 

Pata...aagfig 

The  central  gun  parameters  are  based  on  those  used  in  an  earlier  study 
by  Oberle  [2]  and  are  summarized  in  Tables  2.1  and  2.2.  With  the  exception 
of  Che  plasma  mass,  which  has  been  arbitrarily  set  equal  to  1  g,  and  the 
recognition  of  the  compressibility  of  the  working  fluid,  the  data  of  Table 
2.1  are  as  in  the  study  by  Oberle  [2]  and  correspond  to  the  case  with  water 
as  the  working  fluid.  We  note  the  variation  in  the  effective  energy  as  a 
fraction  of  the  input  electrical  energy  as  the  energy  density  is  varied.  At 
low  energies  and  at  high  energies  the  fraction  of  energy  available  to  do 
work  is  reduced.  For  values  of  the  input  energy  above  the  range  of 
Table  2.2,  a  linear  extrapolation  is  tised  to  get  the  effective  energy. 
Table  2.2  differs  from  the  data  of  Oberle  only  in  chat  a  set  of  values  is 
presented  for  the  case  when  the  input  energy  is  zero.  As  in  the  prevlotis 
study  [15],  we  use  the  Lagrange  pressure  gradient  relations.  The  code 
listing  of  the  nominal  input  data  is  presented  in  Appendix  C  together  with 
Che  complete  solution. 
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In  addition  to  this  nominal  configuration  which  corresponds  to  an  ini¬ 
tial  loading  density  of  46Z,  we  consider  values  of  the  chamber  volume  equal 
to  74.6,  56.0  and  49.8  cm^  corresponding  to  initial  loading  densities  of  60, 
80  and  90X  respectively.  These  allow  us  to  probe  the  influence  of  initial 
ullage  on  the  ballistic  sensitivity  to  the  rate  of  mixing.  In  order  to 
assess  the  Influence  of  the  rate  of  mixing  for  each  loading  density  we 
proceed  as  follows.  First,  we  determine  the  rate  of  plasma  addition  to 
achieve  a  constant  breech  pressure  of  435  MPa,  until  depletion  of  the 
supply,  on  the  assumption  that  the  rate  of  mixing  a(t)  is  proportional  to 
E{(t).  Thus  fli(t)  becomes  equal  to  one  at  the  instant  that  the  supply  of 
plasma  is  depleted.  The  plasma  is  supplied  over  an  Interval  of  approx¬ 
imately  1  ms  at  a  more  or  less  constant  rate.  The  rate  of  delivery  for  each 
of  the  four  loading  densities  is  stimmarized  in  Table  2.3.  The  history  of 
E^(t)  so  determined  is  used  to  perform  runs  in  which  the  rate  of  mixing 
differs  from  the  rate  of  plasma  addition  by  a  constant  factor.  Results  are 
obtained  using  the  thermochemistry  defined  by  Table  2.2  and  with  values 
P  “•  1  and  0  to  assess  the  Influence  of  the  mobility  of  the  unvaporized 
working  fluid.  Results  are  also  obtained  with  the  thermochemistry  frozen  at 
the  10  kJ/g  values  of  Table  2.2  and  with  —  0.  These  results  permit  an 
assessment  of  the  importance  of  variations  in  composition  and  also  provide  a 
basis  for  comparison  with  the  one -dimensional  results  disctissed  in  Chapter 
3.0  since  the  latter  are  performed  subject  to  the  assumptions  of  stationary 
working  fluid  and  constant  mixture  properties. 

Wiim«»r<cal  Results 

Tables  2.4  and  2.5  present  the  influence  of  rate  of  mixing  on  maximum 
breech  pressure  and  muzzle  velocity  for  variable  thermochemistry  and  mobile 
working  fluid.  It  will  be  noted  that  the  maximum  value  of  breech  pressure 
differs  slightly  from  the  target  value  of  435  MPa  %rhen  a  -  for  each  of 
the  initial  loading  densities.  This  is  a  consequence  of  the  tabular 
representation  of  the  ideal  history  of  Ej^  and  could  be  removed  by  using 
greater  resolution  in  the  tabular  ^ta.  However,  for  our  purposes,  it  is 

•  • 

the  variation  with  respect  to  a/E^  which  is  of  interest.  The  same  tabular 
data  are  used  for  all  values  of  a/E^  at  a  fixed  loading  density  and  the 
slight  departure  from  nominal  at  a/E^  —  1  is  not  important. 

We  see  that  the  maximum  pressure  is  relatively  insensitive  to  a  lagging 
of , the  mixing  process.  Maximum  pressure  is  reduced  by  only  about  lOZ  when 
a/E^  -  0.1  for  all  loading  densities.  We  note  that  when  a/E^  -  1,  the  data 
of.  Table  2.2  imply  that  the  effective  energy  will  be  optimal.  When 
a/E^  -  0.1  the  effective  energy  is  a  smaller  fraction  of  the  input  and 
pressurization  is  not  as  intense.  On  the  other  hand,  the  maximum  pressure 
is  much  more  sensitive  to  an  accelerated  rate  of  mixing  with  a/E^  >1.  We 
see  that  in  the  limit  when  the  working  fluid  is  assumed  to  be  entirely 
vaporized  and  mixed  at  the  initial  instant  the  pressure  is  increased 
relative  to  the  ideal  or  proportional  case.  The  low  initial  energy  density 
results  in  low  initial  pressure,  as  is  the  case  for  a/E^  -  0.1.  The  loss  of 
effective  energy  is  more  extreme  at  low  energy  density  than  at  high  energy 
density  as  mfiy  be  seen  from,  the  data  of  Table  2.2.  *^6  difference  between 
the  cases  a/E^  -  0.1  and  a/E^  -  •  stems  from  the  fact  that  in  the  latter  one 
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Table  2.1  Parameters  Used  for  Nominal  Data  Base  (after  Oberle  [2]). 


Chamber  Volume 

97.108  cm’ 

Projectile  Travel 

145  cm 

Bore  Diameter 

14  mm 

Projectile  Mass 

18  g 

Plasma  Energy 

447.8  kJ 

Plasma  Mass 

1  S 

Working  Fluid 

44.78  g  of  H2O 

Density 

1  g/cm’ 

Bulk  Modulus 

5000  MPa 

Derivative  of  Bulk  Modulus 
with  respect  to  pressure 

8. 

Table  2.2 

Thermochemistry  of  Working  Fluid  as  a 
Added  (after  Oberle  [2]). 

Rinctlon  of  Electrical  Energy 

Electrical  Energy  Effective  Energy 

Gamna 

Covolume 

Molecular  Weight 

Eff  En./El.  En. 

Input  (kJ/g) 

(W/g) 

(-) 

(cm’/g) 

(g/mol) 

(-) 

0.0 

0.00 

1.9414 

-1.834 

18.015 

0.000 

3.0 

0.44 

1.9414 

-1.834 

18.015 

0.147 

4.0 

1.63 

1.4000 

-0.342 

18.014 

0.408 

5.0 

2.86 

1.2987 

0.082 

18.014 

0.572 

6.0 

4.08 

1.2557 

0.283 

18.001 

0.680 

7.0 

5.24 

1.2323 

0.404 

17.967 

0.749 

8.0 

6.32 

1.2182 

0.488 

17.896 

0.790 

9.0 

7.31 

1.2093 

0.554 

17.783 

0.812 

10.0 

8.23 

1.2035 

0.609 

17.633 

0.823 

11.0 

9.07 

1.1998 

0.659 

17.450 

0.825 

12.0 

9.85 

1.1975 

0.705 

17.242 

0.821 

13.0 

10.59 

1.1962 

0.748 

17.015 

0.815 

14.0 

11.29 

1.1958 

0.788 

16.775 

0.806 

15.0 

11.95 

1.1959 

0.826 

16.526 

0.797 

16.0 

12.59 

1.1964 

0.863 

16.272 

0.787 

17.0 

13.21 

1.1973 

0.898 

16.015 

0.777 

18.0 

13.81 

1.1985 

0.931 

15.757 

0.767 

19.0 

14.40 

1.1998 

0.964 

15.501 

0.758 

20.0 

14.97 

1.2014 

0.995 

15.248 

0.749 

21.0 

15.55 

1.2030 

1.025 

14.998 

0.740 

22.0 

16.11 

1.2048 

1.054 

14.753 

0.732 

23.0 

16.68 

1.2066 

1.082 

14.512 

0.725 
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Table  2.3  History  of  Plasma  Flux  to  Achieve  Constant  Breech 
Pressure  of  435  MPa 


Fraction  of  Energy  Delivered  to  Mixing  Chamber  (-) 


Time  - 

(msec)  Percent  Initial  Loading  Density  (Z) 


1 

46 

60 

80 

90 

0.0000 

0.2779 

0.1642 

0.0701 

0.0387 

0.0500 

0.2798 

0.1660 

0.0720 

0.0406 

0.1000 

0.2854 

0.1716 

0.0776 

0.0462 

0.1500 

0.2947 

0.1810 

0.0870 

0.0556 

0.2000 

0.3078 

0.1941 

0.1001 

0.0688 

0.2500 

0.3246 

0.2110 

0.1171 

0.0858 

0.3452 

0.2316 

0.1378 

0.1065 

0.3500 

0.3695 

0.2560 

0.1622 

0.1310 

0.3976 

0.2842 

0.1904 

0.1592 

0.4500 

0.4294 

0.3161 

0.2224 

0.1912 

0.4650 

0.3517 

0.2581 

0.2269 

0.5500 

0.5043 

0.3912 

0.2976 

0.2665 

0.5472 

0.4343 

0.3408 

0.3098 

0.6500 

0.5941 

0.4811 

0.3878 

0.3568 

■SVTiTi 

0.6445 

0.5317 

0.4385 

0.4075 

0.7500 

0.6987 

0.5861 

0.4930 

0.4620 

I'H'i'I* 

0.7567 

0.6442 

0.5512 

0.5203 

0.8500  1 

0.8184 

0.7059 

0.6132 

0.5823 

0.8837 

0.7715 

0.6788 

0.6479 

0.9500 

0.9526 

0.8407 

0.7482 

0.7174 

1.0000 

0.9136 

0.8213 

0.7905 

1.0500 

1.0000 

0.9902 

0.8981 

0.8673 

1.0000 

1.0000 

0.9785 

0.9479 

1.1500  I  1.0000  1.0000  1.0000  1.0000 


Table  2.4  Effect  on  Maximum  Breech  Pressure  of  Mixing  Rate  of  Plasma 
with  Working  Fluid.  Thermochemistry  as  in  Table  2.2. 
Unmixed  working  fluid  moves. 


Percent  Pressure  (MPa) 

Initial 
Loading 
Density 
(*) 


Pressure  (MPa) 

P-/Po.i<*) 

Rate  of  Mixing /Rate  of  Plasma  Injection  (-) 
0.1  0.5  1.0  2.0  • 

407.4 

434.3 

435.7 

492.7 

601.9 

1.477 

410.8 

434.9 

436.1 

505.6 

668.5 

1.627 

409.6 

434.7 

437.4 

511.9 

656.4 

1.603 

408.1 

434.4 

438.7 

511.7 

629.3 

1.542 

Table  2.5  Effect  on  Muzzle  Velocity  of  Mixing  Rate  of  Plasma 

with  Working  Fluid.  Thermochemistry  as  In  Table  2.2. 
Unmlxed  working  fluid  moves. 


Percent 

Initial 

Muzzle  Velocity 

(m/s) 

Loaoing 

Density 

Rate 

of  Mixing  /  Rate 

of  Plasma  Injection 

(-) 

(X) 

0.1 

0.5 

1.0 

2.0 

m 

46 

1943.3 

2058.6 

2069.8 

2064 . 7 

2168.6 

60 

1992.9 

2096.6 

2105.2 

2099.6 

2245.4 

80 

2029.2 

2124.2 

2130.1 

2122.1 

2229.1 

90 

2037.8 

2132.1 

2137.4 

2126.8 

2191.7 
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achieves  the  condition  of  high  effective  energy  at  the  time  when  the  plasma 
supply  is  completed.  When  this  occurs  in  the  case  of  instantaneous  mixing 
one  has  the  projectile  travel  reduced  by  comparison  with  the  case  of  propor¬ 
tional  mixing;  the  free  volume  is  reduced  and  the  pressure  then  becomes 
excessive. 

Accordingly,  for  the  case  studied  here,  namely  that  of  water  as  the 
working  fluid,  a  slow  rate  of  mixing  has  only  a  weak  effect  on  ballistics 
stemming  from  an  inefficient  use  of  the  plasma  energy  to  create  propulsive 
gas.  A  rapid  rate  of  mixing,  however,  has  a  stronger  effect  and  is  more 
dangerous  since  it  results  in  increased  pressure.  Of  course,  the  slow  rate 
of  mixing  would  imply  higher  mixture  temperatures  with  the  possibility  of 
greater  heating  of  the  tube.  We  also  note  that  the  sensitivity  to  the  rate 
of  mixing  as  represented  by  the  effect  on  maximum  pressure  is  surprisingly 
indifferent  to  the  initial  loading  density. 

In  Tables  2 . 6  and  2.7  we  probe  the  Influence  of  the  mobility  of  the 
unvaporized  working  fluid.  All  input  data  are  as  for  the  corresponding 
results  in  Tables  2.4  and  2.S  except  that  now  ^  —  0  so  that  the  unvaporized 
working  fluid  is  stationary.  This  implies  that  the  ratios  of  p,  P^^e^ 
p^^  which  were  fixed  in  the  previous  set  of  results  with  »  1,  will  now 
vary  with  a(t).  The  rate  of  plasma  addition  will  be  unchanged.  Initially, 
therefore,  the  spacemean  pressures  will  be  the  same  for  corresponding  cases. 
However,  the  values  of  p^^  will  differ,  the  values  in  Table  2.6  being  lower 
than  those  of  Table  2.4  as  may  be  seen  from  Equations  2.4  and  2.5.  In 
general,  the  trends  with  respect  to  mixing  rate  are  similar  for  the  cases 
p  -  0  and  /J  -  1.  However,  with  ^  -  0  the  sensitivity  is  greater  as 
evidenced  by  the  ratios  of  pressure  for  the  extremes.  Again,  sensitivity  is 
remarkably  indifferent  to  the  value  of  the  loading  density. 

Finally,  in  Tables  2.8  and  2.9  we  examine  the  influence  of  the  varia¬ 
tions  in  thermochemistry  implicit  in  the  preceding  results.  The  results  in 
Tables  2.8  and  2.9  take  the  properties  of  the  mixture  to  be  fixed  at  the 
10  kJ/g  values  of  Table  2.2  and  the  unvaporized  working  fluid  is  assumed  to 
be  stationary,  -  0.  Sensitivity  of  the  maximum  pressure  to  the  rate  of 
mixing  is  reduced  as  evidenced  by  the  ratios  for  the  extreme  cases.  For 
lagging  mixing  rates  there  is  virtually  no  sensitivity  of  the  maximum 
pressure  since  this  is  achieved  at  the  initial  instant.  For  the  accelerated 
mixing  rates  one  has  a  relative  initial  lag  of  projectile  motion  due  to  the 
Increased  pressure  gradient.  As  in  Tables  2.4  and  2.6  this  results  in  a 
lower  free  volume  when  the  plasma  injection  is  complete,  and  hence  higher 
pressure.  Ue  note  in  Tables  2.7  and  2.9  the  relatively  high  velocities 
which  occur  for  the  lagging  mixing  rates.  These  reflect  reduced  pressure 
gradients  when  ^  -  0. 

The  results  of  Table  2.8  are  consistent  with  the  analysis  presented  in 
the  preamble  to  this  section.  Here,  with  constant  thermochemical  proper¬ 
ties,  the  presence  of  initial  ullage,  or  low  loading  density,  acts  to 
desensitize  the  ETG  to  the  variations  in  the  rate  of  mixing.  This  is  the 
case  even  though  p  -  0  and  C/M  -  2.5,  a  relatively  large  ntimber.  Ve  may 
contrast  these  results  with  those  of  Table  2.6,  also  based  on  p  “  0,  where 
we  see  much  greater  sensitivity  to  the  rate  of  mixing  for  all  values  of  the 
initial  loading  density.  Accordingly,  we  may  conclude  that  the  dominant 
influence  in  respect  to  the  sensitivity  exhibited  in  Table  2.6  stems  from 
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Table  2.6  Effect  on  Maximum  Breech  Pressure  of  Mixing  Rate  of  Plasma 
with  Working  Fluid.  Thermochemistry  as  in  Table  2.2. 
Unmixed  working  fluid  stationary. 


P*/P0.1<*> 


2.027 

2.250 

2.181 

2.053 


Percent 

Initial 

Pressure  (MPa) 

Loading 

Density 

(X) 

Rate  of  Mixing /Rate  of  Plasma  Injection  (-) 

0.1  0.5  1.0  2.0  • 

46 

297.0 

365.2 

388.7 

461.8 

601.9 

60 

297.1 

358.7 

376.8 

446.4 

668.5 

80 

301.0 

353.8 

365.7 

429.9 

656.4 

90 

306.5 

353.0 

361.9 

424.5 

629.3 

Table  2 . 7  Effect  on  Muzzle  Velocity  of  Mixing  Rate  of  Plasma 

with  Working  Fluid.  Thermochemistry  as  in  Table  2.2. 
Unmixed  working  fluid  stationary. 


Percent 

Initial 

Muzzle  Velocity 

(m/s) 

Lx>aalng 

Density 

Rate 

of  Mixing  /  Rate 

:  of  Plasma  Injection 

(-) 

(X) 

0.1 

0.5 

1.0 

2.0 

46 

2349.8 

2282.8 

2122.8 

2041.1 

2168.6 

60 

2356.7 

2296.2 

2153.2 

2048 . 6 

22U5.U 

80 

2354.6 

2294.7 

2163.9 

2041 . 1 

2229.1 

90 

2355.3 

2296.8 

2167.3 

2037.0 

2191.7 
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Table  2.8  Effect  on  Maximum  Breech  Pressure  of  Mixing  Rate  of  Plasma 

with  Working  Fluid.  Thermochemistry  fixed  at  10  kJ/g  values 
of  Table  2.2.  Unmixed  working  fluid  stationary. 


Percent 

Initial 

Loading 

Density 

(X) 


Pressure  (MPa) 

Rate  of  Mixing /Rate  of  Plasma  Injection  (-) 

0.1  0.5  1.0  2.0  •• 

382.4 

386.5 

388.1 

412.1 

454.0 

379.8 

378.7 

376.1 

391.3 

472.4 

376.4 

371.2 

364.5 

375.5 

502.0 

373.7 

367.6 

360.1 

371.7 

526.6 

*’-/Po.i<-) 


1.187 

1.244 

1.334 

1.409 


Table  2.9  Effect  on  Muzzle  Velocity  of  Mixing  Rate  of  Plasma  with 
Working  Fluid.  Thermochemistry  fixed  at  10  kJ/g  values 
of  Table  2.2.  Unmixed  working  fluid  stationary. 


Percent  Muzzle  Velocity  (m/s) 

Initial 
Loading 
Density 
(*) 


Muzzle  Velocity 

(m/s) 

Rate 

0.1 

of  Mixing /Rate  of  Plasma  Injection 
0.5  1.0  2.0 

(-) 

m 

2513.7 

2304.0 

2096.6 

2023.6 

2051.0 

2486.1 

2309.3 

2124.7 

2037.2 

2106.4 

2463.9 

2302.1 

2134.4 

2036.3 

2181.0 

2466.2 

2304.2 

2138.2 

2035.7 

2222.3 
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the  dependence  of  the  thermochemical  properties  on  the  Input  energy  density, 
ep.  This  underscores  the  Importance  of  a  proper  representation  of  the 
variation  In  conq>osltlon  of  the  mixture  In  any  future  theoretical  modeling 
of  the  ETG. 

Ue  conclude  this  section  with  the  following  observation.  The  ballistic 
sensitivity  to  the  rate  of  decomposition  exhibited  by  the  present  calcula¬ 
tions  Is  encouraging  from  the  standpoint  of  the  safety  of  the  ETG.  However, 
It  Is  also  evident  that  reliable  modeling  will  require  Information  about  the 
actual  mixing  rates  since  there  Is  a  significant  difference  between 
predictions  based  on  the  two  simplest  hypotheses,  namely  homogeneous  mixing 
and  proportional  mixing. 
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2uJ _ Inverse  Analysis  to  Determine  Rate  of  Deconposltlon  of  Working  Fluid 


As  mentioned  in  Section  2.1,  the  LUMPET  code  contains  an  option  to 
determine  the  rate  of  decomposition  of  the  working  fluid  when  the  pressure 
history  of  the  gun  is  provided  as  an  input.  We  have  also  mentioned  previ¬ 
ously  that  since  the  pressure  is  not  strongly  influenced  by  the  rate  of 
decomposition,  we  must  expect  the  accuracy  of  the  inverse  solution  to  depend 
strongly  on  the  precision  of  the  data  which  characterize  the  plasma  flux, 
the  pressure  history  and  the  thermochemistry. 

To  explore  further  the  influence  of  the  extent  of  decomposition  on  the 
pressure  achieved  for  a  given  supply  of  energy  to  the  working  fluid  in  a 
known  volume,  we  have  written  a  short  code  referred  to  as  PMAP.  We  first 
discuss  this  code  and  present  some  numerical  results  which  illustrate  some 
of  the  difficulties  associated  with  the  inverse  analysis.  Subsequently,  we 
discuss  further  the  inverse  analysis  and  present  some  numerical  results  for 
a  nominal  problem. 

The  PMAP  Code 

We  assume  that  we  are  given  a  total  mass  of  working  fluid  C  in  a  volume 
V.  The  total  energy  available  to  heat  the  vaporized  fraction  of  the  working 
fluid  is  E.  It  follows  from  Equation  2.37  that 


P 


(y  -  l)e^^j(E/Ca) 


(2.43) 


We  neglect  the  initial  ambient  and  the  mass  of  the  plasma  so  that 


_ oC _ 

^  “  V  -  (1  -  a)C/Pp 


(2.44) 


and,  while  we  do  not  consider  the  influence  of  the  energy  of  compression 
stored  in  the  unvaporlzed  working  fluid,  we  do  allow  p  to  be  a  function  of 
pressure  in  accordance  with  Equation  2.23.  We  substitute  2.44  into  2.43  and 
write  the  resulting  equation  for  p  in  the  form 


P 


(r  •  l)e^^^(E/aC) 

1  y  1  -  a  ~ 

a  C  *  ap  ■  ° 

P 


(2.45) 


In  the  previous  section  we  considered  the  case  when  the  thermochemistry  is 
constant  so  that  -  E/aC  and  the  values  of  y  and  b  do  not  depend  on 
E/oC.  Here  we  focus  on  the  case  when  the  thermochemistry  is  variable. 
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The  PMAP  Code  accepts  a  set  of  values  of  E/C  and  C/V,  values  of  and 
Kj  to  support  the  equation  of  state  p_  -  pp(p).  and  a  table  of  values  of 
thenaochemical  data  in  the  form  presented  ^in  Table  2.2.  The  code  then 
determines  p(a)  for  0.025  ^  a  £  1  for  each  pair  of  values  of  E/C  and  C/V. 
Thus  we  can  assess  the  dependence  of  pressure  on  the  fraction  of  the  working 
fluid  which  has  decomposed  for  a  variety  of  values  of  energy  density  and 
loading  density.  A  listing  of  the  code  and  a  detailed  description  of  the 
input  files  are  given  in  Appendix  B. 

Tables  2.10.  2.11  and  2 . 12  present  code  output  for  values  of 
E/C  »  1000,  5000  and  10000  J/g  respectively.  The  working  fluid  is  taken  to 
be  water  with  and  as  in  Table  2,1  and  thermochemical  properties  as  in 
Table  2.2.  The  values  of  C/V  in  each  table  range  from  0.05,  typical  of 
muzzle  exit  conditions,  to  0.5,  typical  of  the  initial  or  early  time 
conditions.  Each  column  presents  the  dependence  of  pressure  on  a  for  a 
fixed  combination  of  E/C  and  C/V. 

An  inspection  of  the  columns  reveals  three  features  which  are  of 
significance  in  respect  to  the  inverse  analysis.  First,  except  for 
E/C  -  10000  J/g  and  C/V  <  0.3,  each  column  shows  at  least  one  maximum  of 
pressure  as  a  function  of  a.  Second,  some  columns,  those  corresponding  to 
the  lower  values  of  E/C  and  the  larger  values  of  C/V,  exhibit  two  maxima. 
Third,  several  coluunns  exhibit  extended  regions  where  the  pressure  depends 
very  weakly  on  a. 

These  observations  have  the  following  implications  for  the  inverse 
analysis.  First,  the  data  used  to  characterize  the  pressure  history  will, 
as  already  anticipated,  have  to  be  very  precise  if  gross  errors  in  the 
determination  of  a  are  to  be  avoided.  Indeed,  it  is  conceivable  that  if  a 
phase  error  were  present  in  respect  to  the  characterizations  of  pressure  and 
plasma  flux,  say  one  lagging  the  other  by  0.1  ms,  conditions  could  easily 
arise  in  which  no  value  of  a  could  be  found  to  match  the  stated  values  of 
energy  and  pressure.  Second,  the  inversion  algorithm  must  deal  with  the 
possibility  of  multiple  roots.  It  is  also  possible  for  those  roots  to  be 
quite  close  together  with  the  result  that  a  bifurcation  occurs  during  the 
time  dependent  inversion  process.  Such  a  condition  will  arise  when  the 
conditions  are  close  to  a  local  maximum  of  pressure  with  respect  to  a. 
Accordingly,  a  principle  of  continuity  may  not  suffice  to  determine  the 
correct  root.  Further  difficulties  in  selecting  the  proper  root  can  stem 
from  experimental  noise  which  may  temporarily  scatter  the  solution  into  an 
inappropriate  branch. 

It  is  emphasized  that  the  problems  associated  with  the  multiplicity  of 
roots  are  a  consequence  of  the  dependence  of  the  thermochemical  properties 
on  the  energy  density.  It  is  obvious  from  an  examination  of  the  equation  of 
state  that  a  unique  solution  for  a  exists  if  the  thermochemical  properties 
are  constant.  However,  the  difficulties  associated  with  the  insensitivity 
of  pressure  to  a  will  remain. 

Inverse  Analysis 

We  consider  two  possible  modes  with  respect  to  the  inverse  analysis. 
In  both  cases  we  assume  that  the  breech  pressure  history  is  given.  In  the 
first  mode  we  assume  that  is  specified  as  a  fixed  value  in  the  input  data. 
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Table  2 . 10  Pressure  as  a  ftinctlon  of  fraction  of  working  fluid  converted 


to  vapor  (a)  and  ratio  of  nass  of  working  fluid  to  chanber 
volume  (C/V) .  Ratio  of  available  energy  to  mass  of  working 
fluid  (E/C)  is  1000  J/g. 


Pressure  (MPa) 

C/V 

(g/cm®) 

« (-) 

.050 

.100 

.200 

.300 

.400 

.500 

.0250 

7.2 

15.1 

34.0 

58.2 

90.0 

133.4 

.0500 

7.9 

16.7 

37.6 

64.3 

99.3 

147.1 

.0750 

8.4 

17.6 

39.5 

67.4 

103.7 

152.6 

.1000 

8.8 

18.5 

41.4 

70.2 

107.6 

157.4 

.1250 

9.0 

19.0 

42.3 

71.6 

109.0 

158.4 

.1500 

9.2 

19.2 

42.6 

71.6 

108.4 

156.2 

.1750 

9.1 

19.1 

42.2 

70.6 

106.0 

151.2 

8.9 

18.6 

40.8 

67.6 

100.6 

141.9 

9.0 

18.8 

40.8 

66.8 

98.0 

136.0 

.2500 

8.4 

17.5 

37.5 

60.9 

88.2 

120.6 

.2750 

10.1 

20.6 

43.4 

68.7 

96.8 

128.0 

9.7 

19.7 

40.6 

62.8 

86.3 

111.4 

.3250 

7.8 

15.7 

31.8 

48.2 

65.0 

82.1 

.3500 

6.9 

13.8 

27.6 

41.5 

55.3 

69.0 

.3750 

6.9 

13.7 

27.3 

40.6 

53.7 

66.7 

.4000 

6.9 

13.6 

26.9 

39.8 

52.3 

64.5 

.4250 

6.8 

13.5 

26.5 

39.0 

50.9 

62.4 

.4500 

6.8 

13.4 

26.2 

38.2 

49.7 

60.5 

.4750 

6.8 

13.3 

25.8 

37.5 

48.4 

58.7 

6.8 

13.3 

25.5 

36.8 

47.3 

57.0 

.5250 

6.7 

13.2 

25.1 

36.1 

46.1 

55.4 

.5500 

6.7 

13.1 

24.8 

35.4 

45.1 

53.9 

.5750 

6.7 

13.0 

24.5 

34.8 

44.1 

52.4 

6.7 

12.9 

24.2 

34.2 

43.1 

51.1 

.6250 

6.6 

12.8 

23.9 

33.6 

42.2 

49.8 

.6500 

6.6 

12.7 

23.6 

33.0 

41.3 

48.5 

.6750 

6.6 

12.7 

23.3 

32.5 

40.4 

47.4 

KZSSl 

6.6 

12.6 

23.1 

32.0 

39.6 

46.2 

.7250 

6.6 

12.5 

22.8 

31.5 

38.8 

45.2 

.7500 

6.5 

12.4 

22.5 

31.0 

38.1 

44.1 

.7750 

6.5 

12.3 

22.3 

30.5 

37.3 

43.2 

6.5 

12.3 

22.0 

30.0 

36.6 

42.2 

.8250 

6.5 

12.2 

21.8 

29.6 

36.0 

41.3 

.8500 

6.4 

12.1 

21.5 

29.1 

35.3 

40.5 

.8750 

6.4 

12.0 

21.3 

28.7 

34.7 

39.7 

6.4 

12.0 

21.1 

28.3 

34.1 

38.9 

.9250 

6.4 

11.9 

20.9 

27.9 

33.5 

38.1 

.9500 

6.4 

11.8 

20.6 

27.5 

32.9 

37.4 

.9750 

6.3 

11.7 

20.4 

27.1 

32.4 

36.7 

1.0000 

6.3 

11.7 

20.2 

26.7 

31.9 

36.0 
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Table  2.11  Pressure  as  a  function  of  fraction  of  working  fluid  converted 
to  vapor  (o)  and  ratio  of  mass  of  working  fluid  to  chamber 
volume  (C/V) .  Ratio  of  available  energy  to  mass  of  working 
fluid  (E/C)  is  5000  J/g. 


« (-) 

Pressure  (MPa) 

.050 

.100 

C/V  (g/cm’) 
.200  .300 

.400 

.500 

,0250 

32.0 

67.4 

150.9 

256.0 

390.4 

565.7 

.0500 

32.9 

69.5 

155.6 

264.0 

584.1 

.0750 

33.9 

71.5 

160.2 

272.0 

415.3 

602.7 

34.9 

73.6 

164.9 

280.1 

427.8 

621.4 

.1250 

35.9 

75.6 

169.6 

288.1 

440.4 

640.2 

36.8 

77.7 

174.3 

296.2 

453.1 

659.3 

.1750 

37.8 

79.8 

179.0 

304.3 

465.8 

678.4 

38.8 

81.9 

183.7 

312.5 

478.6 

697.8 

39.5 

83.3 

186.9 

317.8 

486.7 

709.5 

.2500 

39.7 

83.6 

187.2 

317.7 

485.3 

705.2 

.2750 

40.0 

84.2 

188.2 

318.6 

485.2 

702.6 

40.4 

85.0 

189.4 

319.9 

485.6 

700.6 

.3250 

40.8 

85.8 

190.8 

321.3 

486.2 

698.7 

.3500 

41.2 

86.6 

192,3 

322.9 

487.0 

697.1 

,3750 

41.7 

87.5 

193.6 

324.2 

487.3 

694.8 

.4000 

42.1 

88.2 

194.9 

325.4 

487.4 

692.1 

.4250 

42.5 

89.0 

196.1 

326.4 

487.2 

688.9 

.4500 

42.9 

89.8 

197.2 

327.2 

486.7 

685.3 

,4750 

43.3 

90.4 

198.1 

327.7 

485.6 

680.8 

43.6 

91.0 

198.8 

327.8 

484.0 

675.7 

.5250 

43,9 

91.4 

199.2 

327.5 

481.7 

669.5 

.5500 

44.1 

91.8 

199.4 

326.7 

478.7 

662.5 

.5750 

44.3 

92.1 

199.5 

325.7 

475.5 

655.0 

44.5 

92.3 

199.4 

324.4 

471.6 

646.7 

.6250 

44.6 

92.4 

198.9 

322.5 

467.1 

637.7 

.6500 

44.8 

92.6 

198.6 

320.7 

462.4 

628.2 

.6750 

44.8 

92.6 

197.9 

318.3 

457.0 

617.8 

.7000 

44.8 

92.4 

196.8 

315.4 

450.9 

606.7 

.7250 

44.8 

92.2 

195.8 

312.5 

444.8 

595.4 

.7500 

44.9 

92.2 

194.8 

309.6 

438.4 

583.8 

.7750 

44.8 

91.9 

193.5 

306.1 

431.4 

571.3 

44.7 

91,4 

191.7 

302.0 

423.6 

558.2 

.8250 

44.5 

90.8 

189.6 

297.4 

415.3 

544.6 

.8500 

44.4 

90.6 

188.3 

293.8 

408.0 

531.8 

.8750 

44.5 

90.4 

186.9 

290.0 

400.4 

518.7 

44.4 

90.0 

185.0 

285.6 

392.1 

505.0 

,9250 

44.1 

89.3 

182.7 

280.6 

383.1 

490.6 

.9500 

43.8 

88.4 

180.1 

275.1 

373.6 

475.8 

.9750 

43.4 

87.4 

177.0 

269.1 

363.6 

460.7 

1.0000 

42.9 

86.1 

173.7 

262.7 

353.3 

445.4 
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Table  2.12  Pressure  as  a  function  of  fraction  of  working  fluid  converted 
to  vapor  (a)  and  ratio  of  nass  of  working  fluid  to  chanber 
volume  (C/V) .  Ratio  of  available  energy  to  mass  of  working 
fluid  (E/C)  is  10000  J/g. 


«  (-) 


Pressure  (MPa) 

.050 

.100 

C/V 

.200 

(-) 

.300 

.400 

.500 

62.9 

132.6 

295.8 

498.3 

752.4 

1075.0 

63.9 

134.7 

300.5 

506.6 

765.6 

1095.3 

64.9 

136.7 

305.3 

515.0 

778.9 

1115.8 

65.9 

138.8 

310.1 

523.4 

792.3 

1136.4 

66.8 

140.9 

314.8 

531.8 

805.8 

1157.4 

67.8 

143.0 

319.6 

540.2 

819.4 

1178.5 

68.8 

145.1 

324.4 

548.7 

833.0 

1199.8 

69.8 

147.2 

329.2 

557.2 

846.8 

1221.4 

70.7 

149.3 

334.1 

565.7 

860.7 

1243.3 

71.7 

151.3 

338.9 

574.3 

874.6 

1265.3 

72.7 

153.4 

343.7 

582.9 

888.7 

1287.7 

73.7 

155.5 

348.6 

591.6 

902.8 

1310.3 

74.7 

157.6 

353.4 

600.3 

917.1 

1333.2 

75.6 

159.7 

358.3 

609.0 

931.5 

1356.3 

76.6 

161.8 

363.2 

617.7 

946.0 

1379.7 

77.6 

163.9 

368.1 

626.5 

960.5 

1403.5 

78.6 

166.0 

373.0 

635.4 

975.2 

1427.5 

79.0 

166.8 

374.4 

637.3 

977.4 

1429.2 

79.1 

166.9 

374.1 

635.5 

972.5 

1418.2 

79.3 

167.2 

374.0 

634.0 

967.8 

1407.4 

79.6 

167.6 

374.3 

633.1 

964.0 

1397.8 

79.9 

168.1 

374.7 

632.4 

960.3 

1388.1 

80.2 

168.7 

375.1 

631.8 

956.9 

1378.8 

80.6 

169.3 

375.7 

631.2 

953.5 

1369.6 

80.9 

169.9 

376.3 

630.7 

950.1 

1360.3 

81.3 

170.6 

377.0 

630.3 

946.7 

1350.9 

81.8 

171.3 

377.7 

629.9 

943.5 

1341.7 

82.2 

172.0 

378.4 

629.6 

940.3 

1332.7 

82.6 

172.6 

379.1 

629.1 

936.9 

1323.2 

83.0 

173.3 

379.5 

628.3 

932.9 

1313.2 

83.3 

173.9 

380.0 

627.5 

929.0 

1303.2 

83.7 

174.5 

380.5 

626.5 

924.8 

1292.8 

84.1 

175.1 

380.8 

625.5 

920.6 

1282.4 

84.5 

175.6 

381.2 

624.5 

916.2 

1271.9 

84.8 

176.2 

381.5 

623.3 

911.7 

1261.1 

85.2 

176.7 

381.7 

622.0 

907.1 

1250.4 

85.5 

177.2 

381.8 

620.4 

902.1 

1239.1 

85.8 

177.6 

381.8 

618.7 

896.7 

1227.4 

86.1 

178.0 

381.6 

616.8 

891.3 

1215.7 

86.4 

178.3 

381.4 

614.8 

885.7 

1204.0 
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In  the  second  node  we  do  not  specify  Ins  teed,  the  base  pressure  history 
is  sssuaed  to  be  given.  The  value  of  is  then  peraitted  to  vary  with  tiM 
and  is  deduced  from  the  pressure  gradient  relations  of  Section  2.1.  If  the 
Lagrange  relations  are  used  we  have 


2M 

J* 

H 

1 

_ ^ 

•  Ca  >  a.  -  a 

.  ^base 

i  o 

_ 

C(1  -  a) 


(2.46) 


and,  if  the  chaid>rage  relations  are  used  we  have 


^’br  ’  ^base 

-  aC  -  a.  -  a 

H 

-  a  J  p.  -  a  J  -  b  J 

2  I'^base  It  2  . 

1  o 

I  1 _ 2_± 

C(1  -  a) 


(2.47) 


4^  ilr 

In  Equation  2.47  we  have  a^  •  a./K,  ^  where  a^,  a^ 

and  b  are  defined  by  Equations  2.15  -  2.17. 

Since  it  is  possible  that  Pbase^*^^  contain  experiaental 

noise  we  provide  an  option  to  saooth  the  histories  using  a  high  frequency 
nunerical  filter  developed  by  Shuaan  [16].  The  filter  consists  of  replacing 
the  array  fj  by  where 


-  f j  +  I  (fj^j^  ^J-1  ■  1  <  J  <  N  ,  (2.48) 


and  I,  -  f  1  ,  Ill  ~  f»  where  N  is  the  total  nuod>er  of  data.  As  reconmended  by 
1  1  M  n 

Shuaan,  three  consecutive  passes  are  required  with  w  taking  the  following 
successive  coaplex  values 


wj^  -  0.49965 

(2.49) 

wj  -  -  0.22227  +  0.642401 

» 

(2.50) 

W3  -  -  0.22227  -  0.642401 

• 

(2.51) 

16.  Shuaan,  F.G.  "Numerical  Methods  in  Weather  Predictions: 

II  Smoothing  and  Filtering"  Monthly  Weather  Review  November  1957 
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The  values  of  a(C)  are  determined  from  the  pressure  histories  as 
follows.  Depending  on  the  mode,  p  is  either  given  as  a  fixed  input  or  is 
determined  from  the  values  of  p^^  and  A  trial  value  is  selected  for 
a  and  the  solution  is  advanced  to  the  required  time  level.  The  appropriate 
pressure  gradient  relations  are  used  to  convert  the  input  value  of  breech 
pressure  to  a  corresponding  value  of  spacemean  pressure,  The  computed 
value  of  the  spacemean  pressure,  p(a) ,  is  compared  with  the  value 
corresponding  to  the  input  breech  presstire.  The  value  of  a  is  then  modified 
according  to  a  numerical  representation  of  Newton's  method  (Regula  Falsi)  in 
which  the  derivative  is  computed  numerically  from  the  values  on  successive 
steps .  The  test  function  is  taken  to  be 


f(a)  -  .  1 

Pexp 


(2.52) 


On  the  first  iteration  the  derivative  f(a)  is  not  defined  and  we  simply 
increase  a  by  0.01.  Iteration  continues  until  |f(a)|  is  less  than  some  user 
defined  tolerance,  typically  10"^. 

At  this  point  the  root  so  determined  may  or  may  not  be  accepted.  The 
user  has  the  option  of  specifying  a  nvimber  of  lagging  values  to  which  a 
linear  regression  line  may  be  fitted.  A  total  of  10  values  has  been  found 
appropriate.  If  the  regression  line  is  not  specified  because  the  user  did 
not  request  it  or  because  the  solution  has  not  yet  advanced  sufficiently  in 
time  to  establish  the  required  ntuaber  of  values,  the  first  root  is  accepted 
provided  that  it  is  larger  than  that  at  the  previous  time  step.  If  it  is 
less  than  the  previous  value ,  a  series  of  searches  for  a  larger  root  is 
conducted  with  the  initial  value  being  increased  by  0.01  at  the  beginning  of 
each  search.  The  search  is  terminated  as  soon  as  a  root  is  found  which  is 
larger  than  that  at  the  previous  time. 

In  general  this  approach  is  insufficient  since  noisy  data  may  create 
conditions  in  which  the  solution  for  oi(t)  is  not  monotonlc.  The  requirement 
of  monotonicity  may  push  the  solution  to  an  incorrectly  large  value.  When 
the  regression  line  is  available  a  search  for  a  larger  root  Is  conducted 
according  to  the  preceding  criteria.  However,  the  regression  line  defines 
an  expected  value  and  either  the  smaller  or  the  larger  root  will  be  accepted 
according  to  whichever  is  closer  to  the  expected  value  and  provided  that  the 
trend  of  the  regression  line  is  positive.  If  the  trend  of  the  regression 
line  is  negative,  the  larger  root  is  accepted. 

Ue  note  that  the  search  for  a  second  root  is  only  conducted  if  the 
first  root  is  less  than  the  value  at  the  previous  step,  whether  or  not  the 
regression  line  is  determined.  We  also  tested  an  alternative,  somewhat  more 
complicated,  scheme  in  which  a  search  was  also  made  for  a  smaller  root  if 
the  first  value  of  a  was  larger  than  the  value  at  the  previous  step,  the 
choice  of  the  the  two  roots,  if  indeed  a  second  root  could  be  found,  being 
made  according  to  the  previously  described  criteria.  No  significant 
differences  were  found  for  the  trial  problem  considered  here  and  so  the 
simpler  method  is  to  be  preferred  on  the  basis  of  computational  efficiency. 
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The  trial  value  at  each  step  is  taken  to  be  the  value  at  the  previous 
tlae  step  unless  the  regression  line  is  defined  in  Which  case  we  use  the 
expected  value.  The  slope  and  intercept  of  the  regression  line  are  deter- 
Bined  in  accordance  with  the  standard  fonnilas  [17]. 

Data  Base  and  Numerical  Results 

We  demonstrate  the  inverse  analysis  by  reference  to  a  nominal  problem 
for  which  the  principal  data  are  stumarized  in  Table  2.13.  Ve  see  that  the 
system  parsnwters  are  quite  different  from  those  of  Section  2.2.  The  value 
of  C/M  is  2. 625  and  the  value  of  E/C  -  4762  J/g.  The  themochemical  proper¬ 
ties  of  the  working  fluid  are  described  by  Table  2.2.  The  Lagrange  gradient 
relations  are  used.  The  value  of  ^  is  1.  We  note  that  the  rate  of  decom¬ 
position  of  the  working  fluid  is  exactly  one  half  of  the  rate  of  energy 
flux.  The  complete  input  and  nominal  solution  for  this  problem  are  given  in 
Appendix  D.  It  should  be  noted  that  muzzle  exit  occurs  at  about  3.0  ms  at 
which  time  the  value  of  a  is  0.75.  Thus  the  working  fluid  does  not 
decompose  completely. 

The  computed  histories  of  breech  and  base  pressure  are  used  in  a  varie¬ 
ty  of  ways  to  construct  input  data  to  the  LUMPET  Code  run  in  the  inverse 
mode.  Since  both  pressure  histories  are  available  the  value  of  will  be 
determined  as  part  of  the  inverse  analysis.  All  the  inverse  analyses  use  a 
regression  line  based  on  10  values. 

The  direct  solution  creates  a  table  of  values  of  pressure  for  which  the 
time  interval  is  0.01  ms.  The  pressures  are  printed  to  two  decimal  places 
of  accuracy.  This  complete  table  of  data  is  used  as  input  to  determine  the 
inverse  solution  represented  in  Table  2.14.  We  see  that  the  value  of  fi  is 
determined  almost  exactly.  Even  with  two  decimal  places  of  accxiracy  to 
describe  the  pressure,  however,  the  value  of  a  is  reproduced  with  an  accur¬ 
acy  no  better  than  2.0Z.  We  note  that  the  error  is  largest  at  intermediate 
times  and  in  fact  is  relatively  low  as  the  projectile  nears  the  muzzle. 

In  Table  2.15  we  explore  the  consequences  of  a  sli^tly  less  precise 
representation  of  the  pressure  histories.  We  enter  then  to  one  decimal 
place  of  accuracy.  Since  the  breech  pressure  has  a  maximum  value  of  492  MPa 
and  is  still  equal  to  56  MPa  at  the  muzzle,  this  level  of  accuracy  is  still 
beyond  anything  we  are  liable  to  obtain  experimentally.  We  see  that  small 
errors  now  appear  in  the  values  of  fi,  although  these  are  for  the  most  part 
less  than  IZ.  The  values  of  a  are  <tetermined  with  an  overall  accuracy  which 
is  not  much  worse  than  that  of  the  previous  case,  typically  2  -  3Z.  How¬ 
ever,  large  errors  occur  near  the  muzzle,  of  the  order  of  lOZ.  It  is 
interesting  that  the  inversion  algorithm  stabilizes  itself;  an  error  of  6.0Z 
is  followed  by  one  of  0.9Z  and  an  error  of  10. IZ  is  followed  by  one  of  0.8Z. 
We  also  note  the  non-monotonic  behavior  of  a  towards  the  end  of  the  calcula¬ 
tion. 


17.  Hoel,  P.G.  "Introduce ion  to  M*thematic»l  Statistics" 

John  Wiley  and  Sons  1966 


34 


Table  2.13  ParaneCers  Used  for  Inversion  Study 


Chamber  Volvune 
Projectile  Travel 
Bore  Dlasieter 
Projectile  Mas^ 

Plasma  Energy 
Plasma  Mass 
Working  Fluid  ** 

Density 
Bulk  Modulus 

Derivative  of  Bulk  Modulus 
with  respect  to  Pressure 


444.6  cm^ 
400.0  cm 
4.0  cm 

160  g 

2000000  J 

0.0  g 

420  g  of  H9O 
1  g/cm® 
5000  MPa 
8. 


Plasma  energy  delivered  at  constant  rate  for  2  ms. 
Working  fluid  decomposes  at  constant  rate  for  4  ms. 


Table  2.14  Inverse  Solution  for  a(t)  and  fiit) .  Breech  and  base 

pressure  given  at  every  time  step  to  two  decimal  places 
of  accuracy. 


Time 

(msec) 


Exact 


0.300 

0.350 

0.400 

0.450 

0.500 

0.550 

0.600 

0.650 

0.700 

0.750 


«(-) 

Calculated 


.050 

.101 

.151 

.203 

.253 

0.304 

0.354 

0.406 

0.458 

0.509 

0.557 

0.601 

0.652 

0.699 

0.752 


X  Error 


Exact 


/»(-) 

Calculated 


X  Error 


1.001 
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Table  2.15  Inverse  Solution  for  4i(t)  and  ^(t).  Breech  and  base 

pressure  given  at  every  tine  step  to  one  declnal  place 
of  accuracy. 


Time 

(msec) 


Exact 


0.300 

0.350 

0.400 

0.450 

0.500 

0.550 

0.600 

0.650 

0.700 

0.750 


«(-) 

Calculated 


.051 

.102 

.153 

.205 

.256 

0.307 

0.359 

0.411 

0.464 

0.517 

0.565 

0.636 

0.656 

0.771 

0.756 


1  Error 


2.0 

2.0 

2.0 

2.5 
2.4 

2.3 

2.6 
2.8 
3.1 

3.4 
2.7 
6.0 
0.9 

10.1 

0.8 


Exact 


fii-) 

Calculated 


Z  Error 


1.001 

0.999 

1.001 

0.998 

1.001 

1.003 

0.995 

1.003 

0.994 

1.008 

0.989 


Table  2.16  Inverse  Solution  for  a(t)  and  ^(t).  Breech  and  base 

pressure  given  at  every  time  step  to  zero  decimal  places 
of  accuracy. 


Time 

(msec) 


Exact 


0.300 

0.350 

0.400 

0.450 

0.500 

0.550 

0.600 

0.650 

0.700 

0.750 


a(-) 

Calculated 


.050 
.100 
.149 
.209 
.272 
0.352 
0.427 
0.504 
0.611 
0.684 
0.710 
0.719 
0.948 
0.909 
0.893 


X  Error 


0.0 

0.0 

0.1 

4.5 

8.8 

17.3 

22.0 

26.0 

13.5 

36.8 
29.0 

19.8 

45.8 

29.9 
19.1 


Exact 


1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

l.OOO 


p(-) 

Calculated 


.990 

.003 

.996 

0.995 

0.994 

0.999 

1.008 

0.976 

0.985 

1.035 

1.036 

1.056 

0.731 

0.982 

0.965 


X  Error 
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Ue  carry  the  loss  of  accuracy  one  step  further  In  Table  2.16  in  which 
the  pressures  are  characterized  with  zero  decinal  places.  Now,  even  the 
values  of  become  subject  to  large  sporadic  errors.  At  2.6  ms  the  error  in 
is  26. 9X.  Not  surprisingly,  this  is  accompanied  by  a  very  large  error  in 
a,  namely  45. 8X.  Ue  see  that  the  error  in  a  degrades  rather  quickly  after 
1.0  ms,  being  typically  20X  and  never  less  than  lOX.  Better  results  can  be 
obtained  if  the  pressure  histories  are  smoothed  using  the  algorithm 
described  by  Equations  2.48  -  2.51.  Table  2.17  presents  the  results  of  such 
a  calculation.  We  note  that  the  initial  results  are  less  accurate  than 
those  of  Table  2.16  due  to  the  smearing  of  information.  But  the  later 
results  are  in  general  better,  the  largest  error  in  a  being  23. 2X  while  that 
for  /}  is  5.8X. 

Finally,  in  Table  2.18  we  present  the  results  of  an  inversion  based  on 
a  subset  of  the  complete  pressure  history.  We  only  use  every  fifth  datiim 
from  the  pressure  tables.  Two  decimal  places  of  accuracy  are  retained. 
Values  of  pressure  at  the  missing  time  steps  are  determined  by  linear 
interpolation  of  the  tables.  The  values  of  ft  are  reproduced  faithfully  but 
the  values  of  a  are  in  error  by  10  -  15X  throu^out  most  of  the  calculation. 

In  conclusion,  the  results  of  the  present  study  indicate  that  even  with 
perfect  knowledge  of  the  plasma  flux  and  the  thermochemistry  of  the  working 
fluid,  and  with  the  pressure  specified  to  two  decimal  places,  the  values  of 
a  will  still  be  subject  to  errors  of  the  order  of  2X.  With  regard  to  the 
accuracy  with  which  the  pressure  is  specified  we  should  note  that  the 
effective  limit  is  defined  by  the  accuracy  with  %fhich  the  code  attempts  to 
match  the  data.  In  all  the  present  calculations  this  was  1  part  in  10*  so 
that  specification  of  the  pressure  to  greater  accuracy  would  not  be  expected 
to  have  any  benefit.  If  the  pressure  is  given  with  an  accuracy  of  roughly 
IX,  corresponding  to  the  calculations  with  zero  decimal  places  of  accuracy, 
errors  of  the  order  of  20X  are  to  be  expected. 

The  present  study  has  been  based  on  water  as  the  working  fluid.  It  is 
suggested  that  the  PMAP  Code  be  used  to  screen  possible  alternative  fluids 
to  establish  a  stronger  dependence  of  pressure  on  a.  Even  if  the  candidate 
fluid  were  not  suitable  for  use  in  the  ETC  it  could  nevertheless  be  used  in 
a  test  fixture  to  acquire  valuable  experimental  characterizations  of  the 
rate  of  mixing.  Such  data  would  be  very  useful  as  a  means  of  validating 
multi -dimensional  numerical  simulations  of  the  mixing  process  and,  possibly, 
as  a  means  of  defining  empirical  relations  to  be  used  in  lower  level  codes. 

On  the  other  hand,  working  fluids  which  exhibit  extreme  insensitivity 
of  pressure  to  a  could  be  used  to  determine  the  plasma  flux  from 
measurements  of  pressure  or,  if  the  plasma  flux  were  given,  to  verify  BLAKE 
Code  predictions  of  the  dependence  of  effective  energy  on  electrical  energy 
output . 


Table  2.17  Inverse  Solution  for  a(t}  and  /l(t).  Breech  and  base 

pressure  given  at  every  time  step  to  zero  decimal  places 
of  accuracy.  Breech  and  base  pressure  histories  smoothed. 


Time  a(-) 

(msec)  Exact  Calculated  X  Error 


0.050 

0.100 

0.150 

0.200 

0.250 

0.300 

0.350 

0.400 

0.450 

0.500 

0.550 

0.600 

0.650 

0.700 

0.750 


0.053 

0.105 

0.159 

0.211 

0.266 

0.324 

0.395 

0.473 

0.552 

0.616 

0.656 

0.678 

0.701 

0.806 

0.797 


6.0 

5.0 

6.0 

5.5 

6.4 

8.0 

12.9 

19.5 

22.7 

23.2 

19.3 
13.0 

7.9 

15.1 

6.3 


/*(-) 

Exact  Calculated  X  Error 


.000 

.004 

.998 

.994 

0.994 

0.999 

1.006 

0.984 

0.991 

1.005 

1.021 

1.041 

0.942 

0.982 

0.959 


1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 

1.000 


Table  2.18  Inverse  Solution  for  <s(t)  and  /}(t).  Breech  and  base 

pressure  given  at  every  fifth  time  step  to  two  decimal 
places  of  accuracy. 


Time  oi(-) 

(msec)  Exact  Calculated  Z  Error 


P(-) 

Exact  Calculated  X  Error 


0.250 

0.300 

0.350 

0.400 

0.450 

0.500 

0.550 

0.600 

0.650 

0.700 

0.750 


.056 

.112 

.169 

.224 

0.282 

0.342 

0.402 

0.462 

0.522 

0.582 

0.592 

0.659 

0.668 

0.782 

0.768 


12.0 

12.0 

11.3 
12.0 
12.8 
14.0 
14.8 
15.5 
16.0 

16.4 
7.6 

9.8 

2.8 

11.7 

2.4 


1.000 

1.000 

1.000 


1.000 

1.001 

1.001 


38 


3.0  ONE -DIMENSIONAL  MODEL 


In  the  previous  chapter  we  described  a  lumped  parameter  model  of  the 
ETC  and  used  it  to  demonstrate  weak  sensitivity  of  the  interior  ballistics 
of  the  ETC  to  the  rate  of  decomposition  of  the  working  fluid.  In  the 
present  chapter  we  continue  to  study  the  interaction  between  the  manner  of 
decomposition  of  the  working  fluid  and  the  overall  ballistics  of  the  ETG. 
Our  investigation  is  extended  to  include  the  influence  of  axial  variations 
in  the  state  variables.  Obviously,  a  one -dimensional  continuum  analysis  of 
the  ETG  leaves  much  to  be  desired  since  so  much  of  the  important 
phenomenology,  such  as  the  penetration  of  the  working  fluid  by  the  plasma 
and  the  concomitant  mixing  process,  is  of  an  inherently  two-dimensional 
nature.  Nevertheless,  it  is  of  interest  to  determine  to  what  extent  the 
influence  of  axial  wave  propagation  and  axial  non-uniformity  in  the 
distribution  of  the  working  fluid  and  its  rate  of  decomposition  might 
contribute  to  a  greater  dependency  of  the  overall  ballistics  on  the  details 
of  the  interaction  of  the  plasma  and  the  working  fluid. 

The  role  of  the  one  -  dimens  ional  model  in  the  present  study  is  to  be 
viewed  as  permitting  a  degree  of  simulation  beyond  that  possible  with  the 
lumped  parameter  model.  No  pretense  is  made  of  developing  a  truly 
predictive  model  which,  given  hard  independent  data,  might  be  capable  of 
predicting  maximum  pressure  and  muzzle  velocity  with  any  degree  of  accuracy. 
Accordingly,  we  do  not  develop  a  one -dimens ional  model  from  fundamental 
principles.  Instead,  we  exploit  an  existing  interior  ballistics  code,  XKTC 
[8],  and  make  suitable  minor  modifications  in  order  to  achieve  our  goal.  We 
discuss  the  code  in  Section  3.1.  In  Section  3.2  we  present  some  numerical 
results  for  an  ETG  configuration  similar  to  that  studied  in  Section  2.2. 

lA _ Description  of  Model 

The  XKTC  Code  has  been  developed  over  several  years  and  is  a  general 
purpose  interior  ballistics  code.  It  is  applicable  to  a  wide  variety  of 
solid  propellant  charges  including  both  conventional  and  traveling  charge 
systems.  The  code  supports  a  wide  variety  of  form  functions  of  which  two 
are  of  greatest  relevance  here,  namely  the  monolithic  charge  and  the 
perforated  stick  charge.  Our  discussion  focusses  on  these  particular  charge 
designs.  The  general  details  of  the  XKTC  Code  are  not  reviewed;  the 
interested  reader  is  directed  to  Reference  8  for  a  more  complete  discussion. 
Here  we  simply  discuss  those  aspects  of  XKTC  which  pertain  to  the 
representation  of  the  ETG. 

Figure  3.1  illustrates  the  two  types  of  ETG  configuration  which  are 
admitted  by  the  XKTC  Code  as  revised  in  accordance  with  the  present  task 
objectives.  The  monolithic  charge,  shown  in  Figure  3.1(a),  is  configured 
rather  like  a  solid  propellant  rocket  motor.  It  is  bonded  to  the  sidewall 
of  the  tube  and  has  a  central  port  on  whose  surface  the  combustion  occurs. 
Such  a  configuration  can  be  used  as  an  approximate  representation  of  the 
working  fluid  of  the  ETG,  provided  that  the  working  fluid  is  suitably 
packaged  with  a  central  duct.  In  the  figure  we  show  ullage  inside  the 
monolithic  charge  and  also  between  the  charge  and  the  projectile  base.  We 
refer  to  the  former  as  annular  ullage  and  to  the  latter  as  axial  ullage. 
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Vorkln^  fluid  aa  a  iBonolithlc  charge*  StationaxT'* 
Surface  regreaaion  pTOportional  to  gaa  reloclty* 


(a)  Vorkiug  Fluid  aa  a  Monolithic  Charge. 


Working  fluid  aa  a  perforated  atick 


(b)  Working  fluid  aa  a  Perforated  Stick 


Figure  3*1  XEEC  Bepreaentationa  of  BIG  vith  Annular  Ullage* 
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Evidently  the  annular  ullage  will  increase  in  time  as  the  internal  surface 
of  the  monolithic  charge  decomposes.  The  diameter  of  the  duct,  even  if 
initially  uniform,  may  vary  with  axial  position  in  accordance  with  the  local 
properties  of  the  law  governing  the  rate  of  decomposition.  The  axial  ullage 
will  grow  in  time  due  to  the  motion  of  the  projectile;  however,  decomposi¬ 
tion  of  the  ends  of  the  monolithic  charge  is  not  considered.  Ue  note  that 
although  Figure  3.1(a)  shows  the  rear  face  of  the  monolithic  charge  in 
contact  with  the  breech  face,  the  code  does  admit  the  existence  of  axial 
ullage  to  the  rear  of  the  monolithic  charge.  Such  a  region  would  be  of 
fixed  length  since  the  breechface  is  assumed  to  be  stationary  and  the  rear 
face  of  the  monolithic  charge,  .like  the  forward  face,  is  assumed  not  to 
decompose . 

The  representation  of  Figure  3.1(b)  admits  certain  additional  features. 
The  working  fluid  is  represented  as  a  perforated  stick  which  is  free  to  move 
in  the  axial  direction  in  accordance  with  the  gas -dynamic  forces  exerted  on 
it  by  the  mixture  of  plasma  and  products  of  decomposition.  Annular  ullage 
exists  in  the  region  exterior  to  the  working  fluid  as  well  as  in  the  central 
duct.  Decomposition  can  occur  on  the  outer  and  inner  surfaces,  but  not  at 
the  ends.  The  properties  of  the  gases  in  the  external  region  of  annular 
ullage  are  distinguished  from  those  in  the  interior. 

In  XKTC  the  regions  of  annular  ullage  are  always  represented  according 
to  the  equations  of  motion  for  a  one -dimensional,  unsteady,  compressible 
fluid  with  variations  in  cross-sectional  area  and  with  mass  and  heat 
addition.  The  regions  of  axial  ullage,  at  the  ends  of  the  charge  regions, 
may  be  treated  either  as  continua  or  as  lumped  parameter,  depending  on  their 
lengths.  The  cross-sections  of  the  tube  at  v^ich  the  ends  of  the  charge 
regions  are  located  are  treated  mathematically  as  internal  boundaries, 
permeable  to  the  gas -phase,  but  across  which  all  the  gas -phase  state 
variables  may  Jump  discontinuously  as  a  result  of  the  local  discontinuity  in 
flow  area. 

If  the  working  fluid  is  represented  as  a  monolithic  charge,  the  solid- 
phase  balance  equations  become  trivial  due  to  the  assumption  of  immobility. 
However,  if  the  working  fluid  is  treated  as  a  stick  propellant,  the  motion 
of  the  working  fluid  is  determined  by  the  continuity  and  momentum  equations 
for  an  elastic  rod  of  varying  cross-sectional  area  and  with  mass  transfer  to 
the  ambient. 

The  gas-phase  is  assumed  to  obey  the  covolume  equation  of  state,  as  in 
Chapter  2.0.  However,  the  thermochemical  properties  are  assumed  to  be 
constant;  XKTC  does  not  support  the  BLAKE  Code  formulation  of  variable 
properties . 

The  governing  equations  are  a  system  of  partial  differential  equations 
and  are  solved  using  an  explicit  finite  difference  scheme.  The  discontinu¬ 
ities  at  the  internal  boundaries  defined  by  the  ends  of  the  charge  regions 
are  represented  as  such  in  the  numerical  algorithm. 
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The  representation  of  the  mixing  of  the  plasma  with  the  working  fluid 
is  visualized  as  a  two-step  process.  The  plasma  mixes  with  the  ambient 
gases.  The  motion  of  the  gas,  due  at  first  to  the  non-uniform  heating 
stimulus  by  the  plasma,  and  subsequently  to  the  rarefaction  created  by 
projectile  motion,  induces  decomposition  of  the  working  fluid  throu^  a 
Helmholtz  mechanism. 

The  working  fluid  is  assumed  to  decoiapose  at  a  rate 
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P 


K  lp«  -  pp"pl 
2  (P  +  Pp> 


(3.1) 


where  r  is  the  local  rate  of  surface  regression;  p  is  the  density  of  the 
gas-phase  (the  mixture  of  plasma  and  products  of  decomposition);  p.  is  the 
density  of  the  working  fluid;  u  and  u_  are  respectively  the  velocity  of  the 
gas-phase  and  the  velocity  of  the  working  fluid;  and  is  a  dimensionless 
factor  which  we  refer  to  as  a  wiping  coefficient.  Equation  3.1  has  been 
previously  proposed  in  connection  with  the  regression  of  the  cavity  wall  in 
bulk- loaded  liquid  propellant  guns.  Originally  proposed  to  model 
decomposition  due  to  the  Helmholtz  Instability,  Equation  3.1  can  be  shown  to 
be  closely  related  to  the  Prandtl  theory  of  the  mixing  of  a  tiirbulent  jet 
with  the  ambient  fluid  [18]. 

The  influence  of  the  plasma  is  reflected  in  the  gas -phase  balance 
equations  as  a  source  of  mass,  momentum  and  energy.  As  in  the  lumped 
parameter  model  of  Chapter  2.0,  we  assume  that  the  properties  of  the  plasma 
are  predetermined.  The  plasma  is  assumed  to  be  added  to  the  gas -phase  over 
a  distance  referred  to  as  the  mixing  length.  At  each  time  the  specified 
plasma  flux  is  converted  to  a  distributed  source  term  whose  properties  are 
essentially  uniform  over  the  mixing  length.  To  permit  resolution  by  the 
finite  difference  solver,  the  distribution  is  tapered  linearly  to  zero  over 
the  last  few  mesh  points.  The  magnitude  of  the  source  term  follows  from  the 
value  of  the  flux  according  to  a  simple  quasi- steady  mass  balance. 

Two  representations  of  the  plasma  are  admitted.  In  the  simpler  of  the 
two  we  assume  that  the  plasma  mass  flux  is  given  and  that  the  mixing  length 
is  a  fixed  value  defined  as  an  input  datum.  In  the  second  representation  we 
assume  that  the  state  of  the  plasma  at  the  entrance  to  the  mixing  chamber  is 
completely  specified.  The  mixing  length  is  then  time  dependent  and  follows 
from  the  Prandtl  spreading  rate  for  a  turbulent  jet.  We  take  the  Prandtl 
spreading  rate  for  the  plasma  jet  to  be  [18] 
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18.  Edelwan,  R.B.  "The  Interior  Ballistics  of  Liquid  Propellant 
Guns"  RDA-TR-^408~010 


1974 


42 


where  and  are  the  density  and  velocity  of  the  plasma  jet  as  prescribed 
by  a  table  of  input  data,  and  kj  is  a  dimensionless  coefficient.  The  formal 
similarity  to  Equation  3.1  is*^  obvious.  Ue  assiuie  that  the  Jet  has  an 
initial  diameter  d^ ,  specified  as  input,  and  that  it  is  injected  into  the 
central  ullage  port  whose  diameter  is  d^.  We  then  take  the  mixing  length  to 
be  the  distance  required  for  the  jet  boundary  layer  to  expand  to  the  point 
at  which  it  touches  either  the  wall  of  the  cavity  or  the  centerline  of  the 
jet.  Thus  we  have  the  mixing  length  Ij  defined  by 
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dj.  d^  -  dj] 


(3.3) 


Obviously,  mixing  will  not  be  complete  in  a  distance  I.  which  has  rather  the 
character  of  an  e- folding  length  for  the  mixing  process.  However,  Equation 
3.3  expresses  the  essential  physical  relationship  between  the  jet  parameters 
and  the  rate  of  mixing  with  the  ambient.  The  unknown  coefficient  kj  can  be 
varied  to  compensate  for  the  truncation  inherent  in  Equation  3.3. 

We  note  that  Equation  3.2  is  valid  for  subsonic  jets.  In  fact,  the 
plasma  jet  may  be  underexpanded  and  Equation  3.2  ought  therefore  to  be 
applied  to  conditions  downstream  of  the  shock  system  which  would  be  expected 
near  the  entrance  to  the  mixing  chamber.  Consideration  was  given  to  a 
formulation  similar  to  the  Carfagno  analysis  of  muzzle  blast  as  modified  by 
Hay  and  Einstein  [19].  However,  the  analogy  between  the  two  sets  of 
phenomena  is  not  precise  since  the  muzzle  blast  is  essentially  unconfined. 
In  view  of  the  other  approximations  inherent  in  the  present  application  of 
XKTC,  the  extension  to  model  a  shock  system  in  the  chamber  was  not  thought 
to  be  well  motivated. 

The  goal  of  the  present  study  is  simply  to  assess  the  implications  of 
one -dimensional  axial  structure  of  the  flow  in  respect  to  the  ballistic 
sensitivity  to  the  rate  of  decomposition  of  the  working  fluid.  Accordingly, 
the  results  presented  in  the  next  section  are  based  on  the  simplest  XKTC 
representation  of  the  ETC.  The  working  fluid  is  represented  as  a  monolithic 
charge  and  the  mixing  length  is  taken  to  be  a  fixed  input  datum  for  each 
calculation. 


19.  Hay  I.W.  and  Einstein  ,  S.I.  "Prediction  of  Gun  Muzzle 

Flash"  Proceedings  of  the  lUth  Jannaf  Combustion  Meeting  1977 
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3 ■ 2  Numerical  Results 

The  calculations  presented  here  are  predicated  on  the  same  data  as 
those  in  Section  2.2.  The  mixing  chamber  is  assumed  to  have  a  cylindrical 
section  12.875  cm  in  length  and  a  tapered  section  3  cm  in  length  over  which 
the  diameter  decreases  to  that  of  the  tube.  To  maintain  the  initial  volume 
of  97.1  cm^,  the  diameter  of  the  cylindrical  section  is  set  equal  to  2.944 
cm.  The  working  fluid  is  confined  to  the  cylindrical  section  and  the  cen¬ 
tral  port  has  a  diameter  of  2.059  cm.  To  maintain  consistency  with  Section 
2.2  the  total  plasma  energy  is  set  equal  to  368.6  kJ  which  is  the  effective 
energy  at  10  kJ/g  of  Table  2.2  corresponding  to  an  input  of  447.8  kJ  as  in 
Table  2.1.  The  values  of  y  and  b  are  chosen  consistently  from  Table  2.2. 
The  plasma  is  represented  as  injecting  at  a  uniform  rate  for  0.8  ms.  At  the 
Initial  instant  the  mixing  chamber  is  assumed  to  be  at  atmospheric  pressure 
and  room  temperature. 

As  will  be  seen,  the  solutions  involve  a  great  deal  of  structure  and, 
accordingly,  a  total  of  99  mesh  points  are  used.  XKTC  treats  the  area 
discontinuity  at  the  end  of  the  monolithic  charge  as  an  explicit  internal 
boundary.  Within  each  of  the  two  computational  regions  a  uniform  distribu¬ 
tion  of  mesh  points  is  used,  the  number  being  allocated  to  each  region  in 
proportion  to  length,  except  that  the  monolithic  charge  is  always  allocated 
a  minimum  of  one  half  the  total  number. 

Ue  obtain  results  for  three  values  of  the  plasma  mixing  length,  namely 
2.54  cm,  12.7  cm  and  a  length  which  is  variable  but  always  equal  to  the 
distance  from  the  breech  to  the  base  of  the  projectile.  For  each  mixing 
length  we  obtain  three  solutions  corresponding  to  values  of  k^  -  0.2,  0.4 
and  0.8.  The  results  are  presented  in  Table  3.1  and  in  Figures  3.2  through 
3.13. 


T^Ie  3.1 

Ballistic  Sensitivity  of  EIG  to  FLaaiia  Mixing  Length  and  RatiP  nf  nnmnnpnsitrfnn 
of  Working  Fluid  According  to  One-Dimensional  Gontdium  Analysis 

Plasma  Mixing  Wiping  Coefficient 

a  at  0.4  ms 

a  at  0.8  ms 

Maxinum  Breedi 

Muzzle  Velocity 

Length  (cm] 

(-) 

(-) 

(-) 

FresauacB  (MPa) 

(n/s) 

2.54 

0.2 

0.32 

0.56 

563 

2207 

2.54 

0.4 

0.76 

0.82 

591 

1938 

2.54 

0.8 

0.89 

0.92 

580 

1852 

12.70 

0.2 

0.05 

0.21 

493 

2676 

12.70 

0.4 

0.27 

0.57 

545 

2286 

12.70 

0.8 

0.67 

0.80 

549 

2146 

* 

0.2 

0.01 

0.06 

473 

2906 

* 

0.4 

0.05 

0.29 

468 

2856 

* 

0.8 

0.19 

0.66 

462 

2780 

*  Mixing  length  variable,  taken  to  be  distance  fzva  breedi  to  base  of  projectile 
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In  Table  3.1  we  present  the  maximum  breech  pressure  and  the  nuzzle 
velocity  for  each  of  the  nine  solutions.  We  also  tabulate  the  fraction  of 
the  working  fluid  which  has  been  decomposed  at  0.8  ms,  when  the  plasma 
injection  is  complete,  and  at  0.4  ms,  halfway  through  the  injection  process. 
Naturally,  larger  values  of  the  wiping  coefficient  result  in  larger  values 
of  a  at  a  given  time.  In  general,  the  results  of  Table  3.1  are  consistent 
with  those  of  Section  2.2  based  on  the  lumped  parameter  model.  Higher 
pressures  are  obtained  with  more  rapid  rates  of  decomposition  of  the  working 
fluid. 

Precise  comparisons  between  the  continuum  and  the  lumped  parameter 
results  cannot  be  made  because  the  plasma  injection  r£.te  in  the  latter  were 
chosen  to  give  a  nearly  constant  breech  pressure.  In  the  continuum  results 
the  rate  of  energy  deposition  is  more  gradual  so  that  maximum  pressure 
occurs  at  about  0.4  -  0.8  ms ,  depending  on  the  mixing  parameters .  The 
gradually  rising  pressure  histories  in  the  continuum  solutions  imply  less 
projectile  motion  at  the  time  of  complete  plasma  addition  and  hence  less 
volume  and  more  pressure.  Also,  the  plasma  addition  is  complete  in  0.8  ms 
in  the  continuum  solutions  whereas  in  the  lumped  parameter  solutions  the 
injection  time  was  slightly  longer,  1  ms. 

Figures  3.2,  3.3  and  3.4  present  distributions  of  density,  pressure  and 
velocity  in  the  mixture  of  plasma  and  vaporized  working  fluid  at  various 
times  for  the  shortest  mixing  length  and  the  largest  wiping  coefficient. 
Figure  3 . 5  presents  histories  of  breech  pressure  for  the  shortest  mixing 
length  and  all  three  wiping  coefficients.  An  extraordinary  degree  of 
structure  is  seen  in  the  spacewise  distributions.  The  short  mixing  length 
results  in  a  blast  wave  type  of  response.  A  shock  is  formed  and  undergoes 
reflections  from  the  projectile  base  and  the  face  of  the  monolithic  charge. 
A  highly  structured  flow  is  developed.  The  rate  of  decomposition  is  very 
non-uniform  and  a  slug  of  cool  gas  is  formed  as  the  high  velocity  results  in 
rapid  decomposition  forward  of  the  region  where  heat  addition  due  to  the 
plasma  occurs.  The  wiggles  in  the  pressure  distribution  at  0.2  ms  indicate 
some  strain  on  the  numerical  algorithm.  We  note  that  the  maximum  global 
pressure  is  much  higher  than  the  maximum  breech  pressure  for  this  case.  Yet 
in  spite  of  the  very  non-uniform  distributions  the  pressure  histories  are 
surprisingly  indifferent  to  the  wiping  coefficient  as  shown  in  Figure  3.5. 
Greater  wave  structure  is  seen  with  the  slower  rates  of  decomposition  since 
there  is  less  damping  of  the  initial  explosive  surge. 

Comparable  results  are  presented  for  the  intermediate  plasma  mixing 
length  in  Figures  3.6  through  3.9.  Although  still  highly  structured,  par¬ 
ticularly  as  regards  the  density  distribution,  they  show  less  blast  wave 
structure  than  in  the  previous  case.  The  histories  of  breech  pressure  show 
only  a  modest  sensitivity  to  the  wiping  coefficient,  even  though  the  results 
in  Table  3.1  indicate  marked  differences  in  the  induced  rates  of  decomposi¬ 
tion.  We  see  quite  plainly  the  inflection  due  to  termination  of  the  energy 
supply  at  0.8  ms.  The  lower  pressure  which  occurs  for  the  smallest  value  of 
the  wiping  coefficient  is  associated  with  a  higher  muzzle  velocity,  as  shown 
in  Table  3 . 1  and  reflects  the  reduction  in  the  pressure  gradient  due  to  the 
smaller  mass  of  moving  gas . 


45 


Finally,  in  Figures  3.10  throu^  3.13  we  have  the  conparable  results 
for  the  longest  mixing  length.  The  distributions  are  now  quite  smooth  and 
the  sensitivity  of  the  pressure  history  to  the  rate  of  decomposition  is 
almost  nil,  even  though  the  values  of  a  at  0.8  ms  range  from  0.06  to  0.66  as 
seen  in  Table  3.1  The  density  increases  monotonlcally  through  the  duct  of 
the  monolithic  charge  where  the  working  fluid  is  decomposed  and  added  to  the 
mixture.  In  the  developing  region  of  ullage  behind  the  projectile  base, 
however,  the  density  drops  sharply  as  there  is  no  further  significant  mass 
addition,  only  heating  by  the  plasma. 

Although  a  precise  comparison  with  the  results  of  Section  2.2  is 
precluded  by  the  differences  in  plasma  flux,  it  is  fair  to  say  that  the 
results  of  Table  3.1  are  consistent  with  those  of  Tables  2.8  and  2.9,  for 
which  the  correspondence  is  expected  to  be  greatest.  Even  though  the 
continuum  solutions  reflect  extremely  non>uniform  structure,  they  do  not 
exhibit  ballistic  sensitivity  to  the  rate  of  mixing  in  excess  of  that 
determined  in  the  lumped  parameter  calculations. 
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Figure  3*2  DistrUmtlons  of  density  at  various  tlaes  for 

nixing  length  2.34  on  and  wiping  coefficient  0.8. 
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Figure  3* 


DlBtrllmtloxia  of  pressure  at  Tarious  tiaes  for  plasaa 
airing  length  2.34  ca  and  wiping  ooefficient  0.8. 
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Figure  3*4  l)iatrlbutl<ni8  of  relooity  at  Tarious  times  for  plasma 
mixing  length  2*34  cm  and  vipiag  coefficient  0«8« 
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PBBSSUBS  (MF»  X  10) 


Figure  5*5  Elstories  of  breech  preeenre  for  plasma  airing  length  2*34 
and  three  raluss  of  wiping  ooefficient. 
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Figure  5*6  Dlatributlona  of  (ien8it7  at  Tarloua  times  for  plasma 
mixing  len^fth  12.7  cm  and  wiping  coefficient  0.8* 
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Figure  5*7  Blstrllmtlons  of  pressure  at  yarioua  times  for  plasma 
mixing  length  12*7  cm  and  wiping  coefficient  0.8. 
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3 .8  Distributlona  of  relocity  at  Tarlous  tlmea  for  plasma 
mixing  length  12.7  cm  and  wiping  coefficient  0.8. 
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Figure  5*9  Hlatorlea  of  breeoh  pressure  for  plasna  atxlng  length  12.7  ca 
aiul  rariouB  ralues  of  wiping  coefficient. 
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Figure  3*12  Distrlbutlona  of  relooity  at  rarioua  tloas  for  plaana 
nixing  length  extending  to  projectile  haae  and  wiping 
coefficient  0.8. 
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4.0  TWO-DIMENSIONAL  MODEL 


We  now  consider  Che  formulaclon  of  a  two-dimensional  model  of  Che  ETC. 
Such  a  model  would  permit  an  assessment  of  the  heat  transfer  to  the  tube  as 
well  as  the  overall  ballistic  response.  Section  4.1  contains  a  brief 
summary  of  Che  approach  chat  we  consider  to  be  most  fruitful.  Section  4.2 
presents  the  principal  governing  equations,  a  discussion  of  Che  requirements 
for  closure  of  the  equations,  and  an  appropriate  solution  technique.  In  Che 
present  formulation  two-dimensional  modeling  is  confined  to  the  mixing  cham¬ 
ber  and  the  tube.  The  one -dimensional  model  of  the  capillary  due  to  Powell 
[4]  is  used  to  describe  the  formation  of  Che  plasma.  It  is  the  mixing  of 
the  plasma  with  the  working  fluid  which  is  an  inherently  two-dimensional 
process  and  is  most  naturally  treated  with  a  full  resolution  of  radial  as 
well  as  axial  details  of  Che  flow.  Moreover,  an  important  feature  of  the 
model  is  the  assessment  of  the  heat  transfer  to  the  tube.  This  feature 
absolutely  demands  a  fully  two-dimensional  treatment.  Naturally,  the 
solution  of  the  equations  requires  an  approach  based  on  the  method  of  finite 
differences . 

4.1  Discussion  of  Formulation 

The  flow  in  the  mixing  chamber  is  very  complex  with  as  many  as  three 
distinct  regions  of  highly  structured  flow  each  requiring  detailed  resolu¬ 
tion.  First,  we  must  anticipate  the  formation  of  a  shock  near  the  plasma 
exit  plane  since  the  plasma  is  expected  to  be  underexpanded  at  least  part  of 
Che  time.  The  attendant  flow  will  be  similar  to  chat  associated  with  mvizzle 
blast  and  will  therefore  require  a  fairly  large  number  of  mesh  points  to 
permit  adequate  resolution  by  the  finite  difference  solver.  The  second 
region  of  flow  requiring  careful  resolution  by  the  model  is  that  defined  by 
Che  penetration  of  the  working  fluid  by  Che  plasma  and  the  entrainment  of 
Che  working  fluid  by  the  hot  central  core  flow.  Eventually,  the  central 
cavity  may  penetrate  completely  the  region  occupied  by  the  working  fluid  and 
a  region  of  ullage  may  open  behind  the  projectile.  This  is  of  importance 
since  it  allows  the  heated  gases  easier  access  to  the  tube  wall,  a  problem 
of  great  practical  concern.  The  third  region  requiring  detailed  resolution 
is  Che  boundary  layer  near  the  tube  wall. 

These  requirements  for  resolution  of  flow  details  in  three  regions  — 
near  the  plasma  exit  plane,  in  the  mixing  layer  between  the  core  flow  and 
the  working  fluid,  and  at  the  tube  wall  —  motivate  the  principal  modeling 
restriction,  namely  that  Che  mixture  be  treated  as  locally  homogeneous.  The 
fluid  in  the  mixing  chamber  and  gun  tube  is  viewed  as  a  multicomponent, 
multiphase  mixture  in  which  all  species  and  phases  have  the  same  velocity  at 
each  point  and  all  gas-phase  constituents  have  the  same  temperature.  The 
response  of  the  liquid-phase  is  assumed  to  be  isothermal.  Of  course,  the 
composition  of  the  fluid  is  expected  Co  vary  from  point  Co  point  and  with 
respect  to  time.  These  variations  result  from  convection  of  the  flow, 
turbulent  mixing  and  heat  transfer  and  Che  dependence  of  composition  on 
temperature  and  pressure.  However,  the  assumption  of  local  homogeneity 
implies  that  there  is  no  explicit  internal  boundary  between  the  unvaporized 
working  fluid  and  the  central  region  of  hot  gases.  It  also  means  that  there 
is  no  requirement  Co  formulate  constitutive  laws  governing  Che  formation  of 
droplets  and  their  time  dependent  diameters.  Such  data  are  hard  to  come  by 
and  Che  result  of  including  analysis  to  treat  the  droplet  phase  moving 


independently  from  the  gases  would  be  to  increase  substantially  the  complex¬ 
ity  of  the  code  (i.e.  higher  development  costs  and  longer  run  times)  without 
the  likelihood  of  adding  materially  to  the  predictive  capacity.  It  may  be 
also  pointed  out  that  supercritical  conditions  may  arise  in  which  case 
models  of  the  liquid  phase  as  an  aggregate  of  droplets  are  on  very  tenuous 
ground.  A  useful  discussion  has  been  provided  by  Edelman  [18]  in  which  he 
shows  the  relevance  of  a  turbulent  mixing  law  to  the  analysis  of  the  cavity 
growth  in  the  bulk  loaded  liquid  propellant  gun.  The  algebraic  law  of 
surface  regression  based  on  the  slip  velocity  of  the  liquid  and  gas  phases 
is  shown  to  be  related  logically  to  a  mixing  law  for  a  turbulent  jet  with 
strong  density  gradients.  Finally,  we  note  that  the  local  equilibrium 
approach  has  been  applied  quite  successfully  to  many  types  of  turbulent  jet 
flows  including  that  of  air  penetrating  water  [20,21]. 

The  model  includes  not  only  diffusion  due  to  turbulence  but  also  a 
contribution  due  to  radiative  transport  subject  to  the  simplifying  assump¬ 
tion  that  the  mixture  is  optically  thick.  The  asstimption  of  optical 
thickness  implies  that  the  radiative  transport,  which  is  expected  to  be 
significant,  can  be  expressed  as  a  diffusive  mechanism  [22,23].  A  similar 
approach  has  been  reported  in  studies  of  propulsion  systems  based  on  laser- 
heated  plasmas  [24]  and  by  Powell  in  his  analysis  of  the  plasma  capillary. 

A  two -equation  model  of  turbulence  which  incorporates  differential 
equations  for  the  turbulent  kinetic  energy  and  dissipation  is  thought  to  be 
the  most  logical  choice  in  connection  with  the  two-dimensional  model.  The 
two -equation  model  has  been  used  extensively  in  related  studies  of  both 
turbulent  jets  and  wall  boundary  layers  and  reasonably  good  agreement  with  a 
wide  range  of  conditions  has  been  reported  using  a  fixed  set  of  empirical 
coefficients  [20].  A  one-equation  model,  while  not  as  general,  may  have 
some  relevance  here.  Results  have  been  reported  by  Garloff  and  Helser  in 
connection  with  the  wall  boundary  layer  in  a  gun  tube  [25].  Algebraic 
models  may  be  useful  for  initial  testing.  However,  they  are  not  believed  to 
reflect  accurately  the  conditions  in  a  gun,  even  for  the  simple  case  of  one- 
phase  flow  from  an  initial  all -burnt  condition  in  which  only  the  wall 
boundary  layer  is  present.  They  would  also  necessitate  the  formulation  of 
different  laws  for  the  mixing  region  and  the  wall  boundary  layer. 
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21.  Tross ,  S.R.  "Characteristics  of  a  Submerged  Two-Phase 

Free  Jet"  M.S.  Thesis,  The  Pennsylvania  State  University  197^ 

22.  Oran,  E.S.  and  Boris,  J.P.  "Numerical  Simulation  of 

Reactive  Flow"  Elsevier  1987 

23.  Zel'dovitch,  Ya.B.  and  Raiser,  Yu.P.  "Physics  of  Shock  Waves 

and  High  Temperature  Hydrodynamic  Phenomena"  Academic  Press  1966 

24.  dumb,  R.J.  and  Krier,  H.  "Two-Dimensional  Model  of  Laser- 
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Ve  adopt  the  system  of  equations  presented  by  Garloff  and  Helser  to 
describe  the  motion  of  a  compressible,  viscous,  heat  conducting  and  turbul¬ 
ent  gas  [26].  We  neglect  the  pressure -velocity  correlations.  However,  we 
include  a  governing  equation  for  the  mass  fraction  of  the  fluid  %«hich  is  in 
the  liquid-phase.  It  is  understood  that  all  variables  are  averaged 
quantities  for  the  mixture  of  working  fluid,  its  products  of  decomposition 
and  the  plasma  unless  otherwise  noted. 

Continuity  Equation 


D£ 

Dt 


pV-u 


(4.1) 


where  D/Dt  is  the  convective  derivative. 

Momentum  Equation 

P  5^  -  -  Vp  +  y‘(V  +  ir^)  (4.2) 


Energy  Equation 

^Dt“  ■  *  '^■(Q  +  pt  .  (4.3) 

Turbulent  Kinetic  Energy  Equation 

P  Dt  "  ^  +  0^.  -  pc  -  2p[7(k‘/*))*  (4.4) 


Dissipation  Rate  Equation 

P  k  '  k  ■  2vtV(k‘/*)]*)p  +  7(pjVc)  (4.5) 
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Continuity  of  Liquid  Phase 


DY 


r  -  V(pYV) 


(4.6) 


In  the  foregoing  we  have  used  p,  u,  p,  e,  k  and  c  as  density,  velocity 
vector,  pressure,  internal  energy,  turbulent  kinetic  energy  and  dissipation 
rate  respectively.  We  have  V  and  I|.  as  the  aolecular  and  turbulent  viscous 
stress  tensors,  Q  and  as  the  aolecular  and  turbulent  heat  fluxes.  We 
understand  Q  to  eabed  the  radiative  heat  transport  approximated  as  a 
diffiisive  term.  We  refer  to  Garloff  and  Heiser  [26]  for  further  details  of 
the  turbulence  terms.  We  have  Y  as  the  mass  fraction  in  the  liquid-phase 
and  V  is  the  diffusion  velocity  of  the  liquid.  We  neglect  the  molecular 
contribution  to  this  term  and  resolve  it  as  a  purely  turbulent  phenoswnon. 
The  term  r  in  Equation  4.6  represents  the  rate  of  decosqiosition  of  the 
working  fluid.  We  take  the  gas -phase  to  be  governed  by  a  covolume  equation 
of  state  as  in  Chapter  2.0  with  the  properties  varying  according  to  the 
energy  density.  The  liquid-phase  is  assumed  to  be  isothermal  and  governed 
by  the  eq\iation  of  state  presented  in  Chapter  2.0. 

At  present  a  law  to  describe  r,  the  rate  of  decomposition,  has  not  been 
defined.  Asstming  that  the  rate  of  decomposition  is  a  thermally  dominated 
process,  it  seems  logical  to  propose  that  r  is  proportional  to  the 
difference  between  the  temperatures  of  the  phases.  Possibly,  a  dependence 
on  the  turbulent  velocity  fluctuation  could  be  included  on  the  assumption 
that  the  phases  are  only  in  velocity  equilibrium  in  the  mean.  Further 
attention  to  this  problem  will  be  required  in  future  work. 

An  appropriate  solution  technique  for  the  governing  eqtiations  is 
thought  to  be  a  linearized  alternating  direction  implicit  scheme  of  the  type 
pioneered  by  Beam  and  Warming  [27]  and  by  Briley  and  McDonald  [28]. 
Successful  applications  of  this  techniq\ie  to  the  solution  of  single-phase 
flows  in  gun  tubes  have  been  reported  by  several  authors  [29,30,31]. 
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5.0  CONCLUSIONS 


[1]  A  luaped  paraaeCer  nodel  of  the  ETC  has  been  developed.  The  code 
(LUMPET)  requires  a  characterization  of  the  plasna  flux  and  the 
theraochenical  properties  of  the  aixture  of  plasma  and  working  fluid. 
The  code  can  be  used  to  predict  the  pressure  history  and  projectile 
motion;  to  determine  the  plasma  flux  required  to  achieve  a  constant 
pressure;  to  deduce  the  rate  of  decomposition  of  the  working  fluid 
given  the  pressure  history  in  the  mixing  chafer. 

[2]  Calculations  based  on  a  nominal  data  base  for  the  lumped  parameter  code 
indicate  only  a  weak  dependence  of  the  interior  ballistics  on  the  rate 
of  decomposition  of  the  working  fluid.  Although  the  presence  of  ullage 
is  expected  to  reduce  the  ballistic  sensitivity  to  the  rate  of  decom¬ 
position,  the  variation  of  thermochemical  properties  with  energy 
density  is  found  to  exert  a  dominant  influence. 

[3]  Calculations  based  on  a  nominal  data  base  for  the  Ivunped  parameter  code 
run  in  the  inverse  mode  show  that  very  precise  data  will  be  required 

to  achieve  an  accurate  determination  of  the  rate  of  decomposition  if 
water  is  the  working  fluid.  Even  with  the  values  of  pressure  given  to 
an  accuracy  of  IZ,  errors  of  the  order  of  20Z  can  arise  in  the 
calculated  fraction  of  working  fluid  that  has  been  decomposed. 

[4]  A  code  (PMAP)  has  been  written  to  determine  the  pressure  in  a  chamber 
as  a  function  of  the  fraction  of  the  working  fluid  that  has  decomposed, 
given  valxies  of  the  ratio  of  total  energy  to  total  mass  of  working 
fluid  and  the  ratio  of  total  mass  of  working  fluid  to  total  volume. 

[5]  The  application  of  PMAP  to  the  thermochemical  properties  of  water 
determined  from  the  BLAKE  Code  shows  a  complex  dependence  of  pressure 
on  the  fraction  of  working  fluid  which  has  decomposed.  The  pressure  is 
non-monotonic  tuider  some  conditions;  one  or  more  maxima  may  be 
observed.  Accordingly,  the  inverse  analysis  in  which  the  rate  of 
decomposition  is  to  be  determined  from  the  pressure  history  may 
encounter  multiple  roots. 

[6]  An  existing  one -dimensional,  two-phase  flow  model  of  solid  propellant 
interior  ballistics  (XKTC)  has  been  modified  to  permit  simulations  of  a 
class  of  ETC  configurations  in  %fhich  a  central  ullage  port  is  present. 
As  with  the  lumped  parameter  models,  it  is  assumed  that  the  plasma  flux 
is  provided  as  input  data. 

[7]  The  application  of  XKTC  to  a  nominal  data  base  revealed  no  ballistic 
sensitivity  to  the  rate  of  decomposition  of  the  working  fluid  in  excess 
of  that  determined  using  the  lumped  parameter  model.  This  was  true  in 
spite  of  the  fact  that  the  contlnuiim  solutions  contained  an  extra¬ 
ordinary  level  of  structure.  However,  the  non-uniformity  was  such  as 
to  suggest  that  continuum  modeling  will  be  necessary  to  provide 
reliable  estimates  of  the  pressure  gradient. 

[8]  A  two-dimensional  model  of  the  ETC  has  been  formulated. 
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NOMENCLATURE 


A 

A(z) 

b 

C 


Ea 

Ei(C) 

®^kin 


Pre- exponent  in  Arrhenius  reaction  rate  law 

Cross-sectional  area  at  distance  z  froa  rear  of  mixing  chasd>er 

Bore  area 

Covoluae 

Total  mass  of  working  fluid 

Heat  capacities  at  constant  pressure  and  volume 

Turbulence  parameters 

Diameter  of  cavity  in  working  fluid 

Diameter  of  plasma  Jet  at  capillary  exit  plane 

Activation  energy 

Total  plasma  energy  delivered  to  mixing  chaad>er  at  time  t 

Kinetic  energy  of  working  fluid  and  mixture  of  products  of 
decomposition  and  plasma 

Internal  energy  of  initial  ambient  in  mixing  chamber 

Internal  energy 

Effective  internal  energy 

Energy  of  compression  of  working  fluid 

Turbulence  parameters 

Constant  to  reconcile  units  of  measurement 
Heat  of  formation 

Integrals  associated  with  chambrage  analysis  of 
pressure  gradient 

Turbulent  kinetic  energy 

Dimensionaless  coefficient  to  control  decomposition  of 
working  fluid  according  to  Helmholtz  mechanism  in  XKTC  Code 

Dimensionless  coefficient  to  contol  rate  of  mixing  of  plasma  Jet 
with  ambient  in  XKTC  Code 

Bulk  modulus  of  working  fluid 
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Mixing  length  of  plasma  Jet 
Mass  of  projectile 
Molecular  weight 

Total  mass  of  plasma  delivered  to  mixing  chamber  at  tliM  t 

Mass  of  initial  aid>ient  in  mixing  chamber 

Pressure 

Molecular  heat  flux  vector 

Turbulent  heat  flux  vector 

Heat  loss  to  chamber  and  tube 

Universal  gas  constant 

Rate  of  decomposition  of  working  fluid 

Rate  of  mixing  of  plasma  jet 

Temperature 

Time 

Velocity  of  gas 
Velocity  of  working  fluid 

Velocity  of  plasma  jet  at  capillary  exit  plane 
In  Chapter  2.0,  volume  of  gas 

In  Chapter  4.0,  diffusion  velocity  of  working  fluid 
Volume  of  mixing  chamber 

Volume  contained  between  rear  face  of  mixing  chamber  and 
position  z 

Projectile  velocity 

Work  done  against  resistive  forces  acting  on  projectile 
Displacement  of  projectile 

In  Chapter  2.0,  mass  fraction  of  intermediate  species  of 
decomposition 

In  Chapter  4.0,  mass  fraction  of  working  fluid 
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Greek  Symbols 

Fraction  of  working  fluid  that  has  decomposed  at  time  t 
Lagrange  coupling  coefficient 
lUitio  of  specific  heats 

Dissipation  rate  of  turbulent  kinetic  energy 
Viscosity 

Effective  viscosity  in  k-eqxiation 

Effective  viscosity  in  c>eqtiation 

In  Chapter  2.0,  order  of  reaction 
In  Chapter  4.0,  kinematic  viscosity 

Molecular  viscous  stress  tensor 

Turbulent  viscous  stress  tensor 

Density  of  gas 

Density  of  plasma  Jet  at  capillary  exit  plane 
Density  of  working  fluid 

In  Chapter  2.0,  parameter  set  equal  to  zero  or  one  depending  on 
representation  of  thermochemistry  of  products  of  decomposition  of 
working  fluid 

In  Chapter  4.0,  molecular  mechanical  dissipation  function 
Turbulent  mechanical  dissipation  function 


INTEN” r ''NAL.  y  -.EFT  BLANK. 


72 


APPENDIX  A 

LUMPET:  Code  Listing  and  Description  of  Input 
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Intent lONALLy  left  blank. 


C*MFST  $ST0RAGE:2 

C  PROGRAM  TO  PERFORM  LUMPED  PARAMETER  SIMULATION  OF  ELECTROTHERMAL 

C  GUN  INTERIOR  BALLISTICS. 

C 

C  VERSION  1.0  (DECEMBER  1989)  WRITTEN  FOR  F77  CCXIPILER  ON  CRAY. 

C 

C  **************************  it  a****************************** ******** 

C  CODE  MAY  BE  CONVERTED  TO  VERSION  COMPATIBLE  WITH  MICROSOFT  F77 

C  CCmPILER  BY  MEANS  OF  GLOBAL  TEXT  EDITING  CCNOfAND  OF  THE  FORM: 

C 

C  REPLACE  'C*MSFT  '  WITH  '  ». 

C 

C  ***************************************'*****'**********'*******"**■*** 

C  MICROSOFT  F77  VERSION  PROMPTS  AT  TERMINAL  FOR  INPUT  DATA  FILE  NAME 

C  (FILNAM)  AND  OUTPUT  FILE  NAME  (OUTFIL) .  THE  OUTPUT  WILL  BE 

C  DIRECTED  TO  THE  PRINTER  IF  OUTFIL  IS  ENTERED  AS  PRN. 

C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


DESCRIPTION  OF  INPUT  DATA  (ASSUMED  TO  BE  LOCATED  IN  FILE  LABELED 

FILNAM) : 

FILE  1.0:  "TITLE"  (A80)  ONE  CARD. 

TITLE  -  PROBLEM  TITLE.  UP  TO  80  CHARACTERS. 

FILE  2.0:  "INTEGRATION  PARAMETERS"  (6F10.0)  ONE  CARD. 

DTNOM  -  INTEGRATION  TIME  STEP  (MSEC). 

PTOL  •  CONVERGENCE  TOLERANCE  IN  SOLUTION  FOR  PRESSURE  (-)• 

DEFAULT  IS  0.0001.  VALUES  SMALLER  THAN  THIS  MAY 
NOT  BE  ADVISABLE  ON  A  PC  UNLESS  THE  PROGRAM  IS 
CONVERTED  TO  DOUBLE  PRECISION.  FAILURE  TO  CONVERGE 
TO  WITHIN  THE  VALUE  OF  PTOL  RESULTS  IN  TERMINATION  OF 
THE  CALCULATION. 

PTOLC  •  CONVERGENCE  TOLERANCE  IN  DETERMINATION  OF  CONSTANT 
PRESSURE  SOLUTION  (-)•  DEFAULT  IS  0.0001.  ONLY 
REQUIRED  IF  A  CONSTANT  PRESSURE  SOLUTION  IS  BEING 
SOUGHT.  FAILURE  TO  CONVERGE  TO  WITHIN  THE  VALUE  OF 
PTOLC  RESULTS  IN  TERMINATION  OF  THE  CALCULATION. 

TTOL  •  CONVERGENCE  TOLERANCE  IN  SOLUTION  FOR  TEMPERATURE  (-)• 

DEFAULT  IS  0.0001.  ONLY  REQUIRED  IF  FINITE  RATE 
CHEMISTRY  IS  CONSIDERED.  FAILURE  TO  CONVERGE 
TO  WITHIN  THE  VALUE  OF  TTOL  RESULTS  IN  TERMINATION 
OF  THE  CALCULATION. 

ATOLl  -  FIRST  CONVERGENCE  TOLERANCE  IN  SOLUTION  FOR  FRACTION 
OF  WORKING  FLUID  WHICH  HAS  DECOMPOSED.  DEFAULT  IS 
0.0001.  ONLY  REQUIRED  IF  MATCHING  TO  OBSERVED 
BALLISTICS  IS  USED  AS  A  MEANS  OF  DETERMINING  THE 
RATE  OF  DECOMPOSITION.  FAILURE  TO  CONVERGE 
TO  WITHIN  THE  VALUE  OF  ATOLl  RESULTS  IN  CONTINUATION 
OF  THE  CALCULATION  AND  PRINTING  OF  AN  ERROR  MESSAGE. 
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AT0L2  -  SECOND  CONVERGENCE  TOLERANCE  IN  SOLUTKMl  FOR  FRACTKM 
OF  tfORKING  FLUID  WHICH  HAS  DECOMPOSED.  DEFAULT  IS 
0.01.  ONLY  REQUIRED  IF  MATCHING  TO  OBSERVED 
BALLISTICS  IS  USED  AS  A  MEANS  OF  DETERMINING  THE  RATE 
OF  DECOMPOSITION.  FAILURE  TO  CONVERGE  TO  WITHIN  THE 
VALUE  OF  ATOL2  RESULTS  IN  TERMINATION  OF  THE 
CALCULATION. 

FILE  3.0:  "PRINT  CYCLE  COUNTER"  (15)  ONE  CARD. 

NPRT  -  NUMBER  OF  STEPS  BETWEEN  SUCCESSIVE  TABULATIONS  OF 

SOLUTION. 

FILE  4.0:  "GENERAL  DATA"  (8F10.0)  ONE  CARD. 

PRM  -  PROJECTILE  MASS  (G) . 

PRTRV  -  PROJECTILE  TRAVEL  (CM). 

DB  -  BORE  DIAMETER  (CM). 

VO  -  VOLUME  OF  MIXING  CHAMBER  (CC) . 

TEMPO  -*  INITIAL  TEMPERATURE  OF  GAS  IN  MIXING  CHAMBER  (K)  . 

PO  -  INITIAL  PRESSURE  OF  GAS  IN  MIXING  CHAMBER  (MPA) . 

PBSK  -  DESIRED  VALUE  OF  BREECH  PRESSURE  FOR  CONSTANT 

PRESSURE  BALLISTICS  (MPA).  IF  PBSK  IS  ENTERED 
AS  ZERO,  A  CONSTANT  PRESSURE  SOLUTION  IS  NOT  IMPOSED. 
THE  PRESSURE  FOLLOWS  FROM  THE  RATE  OF  DEPOSITION  OF 
PLASMA  ENERGY  SPECIFIED  IN  FILE  10.0.  IF  PBSK  IS 
ENTERED  AS  NON-ZERO,  THE  RATE  OF  PLASMA  ADDITION  TO 
ACHIEVE  THE  DESIRED  CONSTANT  PRESSURE  IS  DETERMINED 
AS  PART  OF  THE  SOLUTION,  SUBJECT  TO  THE  (XRISTRAINT 
THAT  THE  TOTAL  PLASMA  ENERGY  BE  LIMITED  TO  THE  VALUE 
TOTPE  DEFINED  IN  FILE  9.0.  IT  IS  ALSO  ASSUMED  IN  THIS 
CASE  THAT  THE  RATE  OF  DECOMPOSITION  OF  THE  WORKING 
FLUID  IS  PROPORTIONAL  TO  THE  RATE  OF  ADDITION  OF 
PLASMA  ENERGY  IN  ACCORDANCE  WITH  THE  VALUE  OF  RATE. 

RATE  •  FRACTIONAL  RATE  OF  DECOfPOSITION  OF  WORKING  FLUID 

DIVIDED  BY  FRACTIONAL  RATE  OF  DEPOSITION  OF  PLASMA 
ENERGY  (•).  IF  RATE  IS  ENTERED  AS  ZERO,  THE  RATE  OF 
DECOMPOSITION  OF  THE  WORKING  FLUID  IS  ASSUMED  TO  BE 
DEFINED  BY  THE  DATA  OF  TABLE  12.0.  IF  RATE  IS  NOT 
ENTERED  AS  ZERO,  THE  RATE  OF  DECOMPOSITION  IS  ASSUMED 
TO  BE  PROPORTIONAL  TO  THE  FRACTIONAL  RATE  OF 
DEPOSITION  OF  PLASMA  ENERGY  AS  DEFINED  IN  FILE  10.0 
SUBJECT  TO  THE  CONSTRAINT  THAT  THE  TOTAL  FRACTION  OF 
THE  WORKING  FLUID  WHICH  HAS  DECOMPOSED  AT  ANY  TIME 
BE  LESS  THAN  OR  EQUAL  TO  ONE.  A  VALUE  OF  RATE  -  I 
PRODUCES  PROPORTIONAL  DECOMPOSITION.  A  VALUE  >  1 
PRODUCES  ACCELERATED  DECOMPOSITION  AND  A  VALUE  <  1 
PRODUCES  LAGGING  DECOMPOSITION. 

HOWEVER,  IF  THE  DATA  OF  FILE  16  ARE  PROVIDED,  AN 
INVERSE  ANALYSIS  IS  PERFORMED  IN  WHICH  THE  RATE  OF 
DECOMPOSITION  IS  DETERMINED  TO  MATCH  THE  OBSERVED 
PRESSURE  HISTORY.  IN  THIS  CASE  THE  VALUE  OF  RATE 
IS  IGNORED. 
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FILE  4.1: 
NSTA 


NSPLT 


FILE  4.21 
N.B.  THIS 
ZA(1) 
RA<1) 
ZA(2) 


RA(NSTA) 

NOTES : 

FILE  5.0: 
NRES 


NSHK 


THE  INITIAL  CONDITION  WILL  (XILY  YIELD  TEMPO  AMD  PO 
IF  THE  PLASMA  IS  SPECIFIED  AS  HAVING  ZERO  TRANSFER 
AT  THE  INITIAL  TIME.  TEMPO  AND  PO  MAY  BE  THOUGHT  OF 
AS  SPECIFYING  THE  AMBIENT  (XMIDITIOM  JUST  BEFORE  THE 
TRANSFER  OF  PLASMA  BEGINS. 

"NUMBER  OF  DATA  TO  DESCRIBE  CHAMBER  GEOMETRY” 

(215)  ONE  CARD. 

-  NUMBER  OF  PAIRS  OF  DATA  IN  FILE  4.2.  IF  NSTA  IS 
ZERO.  FILE  4.2  IS  NOT  REQUIRED. 

NSTA  ALSO  CONTROLS  THE  MODEL  OF  THE  PRESSURE  GRADIENT. 
IF  NSTA  IS  ZERO  THE  STANDARD  LAGRANGE  RELATIONS  ARE 
USED  TO  RELATE  THE  BREECH.  BASE  AND  SPACEMEAM 
PRESSURES  AMD  TO  DETERMINE  THE  KINETIC  ENERGY  OF 
THE  MOVING  FLUID. 

IF  NSTA  IS  GREATER  THAN  OR  EQUAL  TO  TWO.  THE 
CHAMBRAGE  RELATIONS  AS  DESCRIBED  BY  ROBBINS  ET  AL 
ARE  USED. 

NSTA  MUST  NOT  EXCEED  TWENTY. 

-  NUMBER  OF  INCREMENTS  FOR  DETERMINATION  OF  CHAMBER 
INTEGRALS  IN  CHAMBRAGE  ANALYSIS.  EACH  INTERVAL 
(ZA(I-l).ZA(I))  IS  SPLIT  INTO  NSPLT  INCREMENTS. 

DEFAULT  IS  ONE. 

"CHAMBER  GEOMETRY"  (2F10.0)  NSTA  CARDS. 

FILE  IS  REQUIRED  IF  AMD  ONLY  IF  NSTA  >  1. 

-  FIRST  DISTANCE  FRCRf  REAR  FACE  OF  CHAMBER  (CM) . 

-  CORRESPONDING  RADIUS  OF  CHAMBER  (CM). 

-  SECOND  DISTANCE.  STARTS  A  NEW  CARD. 


THESE  DATA  ARE  LINEARLY  INTERPOLATED. 

"RESISTANCE  PARAMETERS”  (215)  ONE  CARD. 

-  NUMBER  OF  PAIRS  OF  DATA  TO  CHARACTERIZE  BORE 
RESISTANCE.  IF  NRES  IS  ZERO.  THE  BORE  RESISTANCE 
IS  ASSUMED  TO  BE  ZERO.  IF  NRES  IS  GREATER  THAN 
ZERO.  FILE  6.0  IS  REQUIRED.  NRES  NAY  NOT  EXCEED 
20. 

-  IF  NSHK  IS  ZERO,  THE  RESISTANCE  DUE  TO  THE  AIR  IN 
FRONT  OF  THE  PROJECTILE  IS  NOT  CONSIDERED. 

IF  NSHK  IS  ONE,  THE  PRESSURE  OF  THE  AIR  IN  FRONT 
OF  THE  PROJECTILE  IS  CONFUTED  USING  STEADY  STATE 
SHOCKWAVE  THEORY.  THE  AIR  PRESSURE  IS  ADDED  TO 
THE  BORE  RESISTANCE  PRESSURE. 
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FILE  6.0;  "  BORE  RESISTANCE  AS  A  FUNCTION  OF  PROJECTILE  TRAVEL" 
(2F10.0)  NRES  CARDS. 

N.B.  THIS  FILE  IS  REQUIRED  IF  AND  ONLY  IF  NRES  >  0. 

RESX(l)  -  FIRST  VALUE  OF  PROJECTILE  DISPLACEMENT  (CM). 

RESP(l)  -  CORRESPONDING  VALUE  OF  RESISTIVE  PRESSURE  (MPA). 
RESX(2)  •  SECOND  VALUE  OF  DISPLACEMENT.  STARTS  A  NEW  CARD. 


RESP(NRES) 

NOTES:  OUTSIDE  THE  TABLE  RANGE  THE  FIRST  OR  LAST  VALUE  IS 

USED.  THUS  A  CONSTANT  BORE  RESISTANCE  CAN  BE 
CHARACTERIZED  WITH  ONE  PAIR  OF  DATA. 

INSIDE  THE  TABLE  RANGE  THE  DATA  ARE  LINEARLY 
INTERPOLATED. 

FILE  7.0:  -PROPERTIES  OF  UNDECOMPOSED  WORKING  FLUID" 

(6F10.0)  ONE  CARD. 

C  -  TOTAL  MASS  OF  WORKING  FLUID  (G) . 

RHOPO  -  DENSITY  OF  WORKING  FLUID  AT  ONE  ATMOSPHERE  (G/CC) . 

AKl  •  BULK  MODULUS  AT  ONE  ATMOSPHERE  (MPA) . 

AK2  -  DERIVATIVE  OF  MODULUS  WITH  RESPECT  TO  PRESSURE  (-)• 

QV  -  HEAT  OF  VAPORIZATION  OF  WORKING  FLUID  (J/G) . 

BET  -  LAGRANGE  COUPLING  COEFFICIENT  (-).  MAY  BE  ENTERED  AS 

ANY  VALUE  BETWEEN  ZERO  AND  ONE.  IF  BET  IS  ZERO,  THE 
UNDECOMPOSED  WORKING  FLUID  IS  ASSUMED  TO  BE  AT  REST 
AT  ALL  TIMES.  IF  BET  IS  ONE,  THE  UNDECOMPOSED  FLUID 
IS  ASSUMED  TO  HAVE  THE  SAME  VELOCITY  PROFILE 
(LAGRANGE)  AS  THE  MIXTURE  OF  PLASMA  AND  DECOMPOSED 
FLUID. 

NOTES:  QV  SHOULD  BE  SET  EQUAL  TO  ZERO  IF  THE  PROPERTIES 

OF  THE  DECOMPOSED  WORKING  FLUID  ARE  CHARACTERIZED  BY 
THE  PARAMETRIC  BLAKE  CODE  DATA  OF  FILE  14.0.  THE  DATA 
OF  FILE  14.0  AUTOMATICALLY  INCLUDE  THE  ENERGY  REQUIRED 
TO  DECOMPOSE  THE  WORKING  FLUID. 

FILE  7.1:  "PROPERTIES  OF  FINAL  DECOMPOSITION  PRODUCTS" 

(3F10.0)  ONE  CARD. 

B  -  COVOLUME  OF  FINAL  DECOMPOSITION  PRODUCTS  (CC/G) . 

GAM  -  RATIO  OF  SPECIFIC  HEATS  OF  FINAL  DECOMPOSITION 

PRODUCTS  (-). 

GMOL  -  MOLECULAR  WEIGHT  OF  FINAL  DECCMtPOSITION 

PRODUCTS  (G/GMOL). 

NOTES:  IF  THE  DATA  OF  FILE  14.0  ARE  PROVIDED  THE  VALUES  OF 

B,  GAM  AND  BET  ARE  ONLY  USED  TO  DETERMINE  THE  INITIAL 
CONDITIONS . 
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FILE  7.2:  "PROPERTIES  OF  INTERMEDIATE  DECOMPOSITION  PRODUCTS " 

(E10.2.7F10.0)  ONE  CARD. 

ARCB  -  PRE- EXPONENTIAL  FACTOR  IN  ARRHENIUS  RATE  LAW  (UNITS 

YIELD  G/CC-SEC).  IF  ARCB-O,  IT  IS  ASSUMED  THAT 
IITERMEDIATE  PRODUCTS  ARE  NOT  FORMED  AND  THE 
REMAINING  DATA  IN  THIS  FILE  ARE  NOT  REQUIRED. 

ARXB  -  TEMPERATURE  EXPONENT  IN  ARRHENIUS  LAW  (•). 

AREB  -  ACTIVATION  ENERGY  IN  ARRHENIUS  RATE  LAW  (J/GMOL)  . 

AROB  -  ORDER  OF  REACTION  WITH  RESPECT  TO  CONCENTRATION  OF 

INTERMEDIATE  PRODUCTS  (-). 

ECHB  -  CHEMICAL  ENERGY  RELEASED  BY  TRANSFORMATION  OF 

INITIAL  PRODUCTS  OF  DECOMPOSITION  INTO  FINAL 
EQUILIBRIUM  PRODUCTS  (J/G) .  ECHB  IS  POSITIVE  IF 
THE  REACTION  TO  FINAL  PRODUCTS  IS  EXOTHERMIC. 

BI  -  COVOLUME  OF  INTERMEDIATE  DECOMPOSITION  PRODUCTS  (CC/G) 

GAMI  -  RATIO  OF  SPECIFIC  HEATS  OF  INTERMEDIATE  DECOMPOSITION 

PRODUCTS  (-). 

GMOLI  -  MOLECULAR  WEIGHT  OF  INTERMEDIATE  DECOMPOSITION 
PRODUCTS  (G/GMOL). 

NOTES:  ECHB  IS  NOT  USED  IF  THE  PROPERTIES  OF  THE  FINAL 

PRODUCTS  ARE  DEFINED  BY  THE  DATA  OF  FILE  14.0. 

IT  IS  ASSUMED  IN  THAT  CASE  THAT  THE  HEAT  OF 
REACTION  IS  INCORPORATED  INTO  THE  EFFECTIVE 
INTERNAL  ENERGY  IN  THE  SAME  HAY  THAT  THE  HEAT 
OF  VAPORIZATION  IS  CAPTURED  BY  THE  EQUILIBRIUM 
ANALYSIS  OF  THE  BLAKE  CODE.  FILE  14.0  PROVIDES 
A  TABULATION  OF  THE  DIFFERENCE  IN  THE  HEATS  OF 
FORMATION  OF  THE  INTERMEDIATE  AND  FINAL  PRODUCTS 
WHICH  SUPERSEDE  THE  VALUE  OF  ECHB. 

IT  IS  ASSUMED  THAT  THE  RATE  OF  DECOMPOSITION  OF 
THE  WORKING  FLUID  IS  GIVEN  IN  TABULAR  FORM  WHEN 
FINITE  RATE  CHEMISTRY  IS  CONSIDERED. 

FILE  8.0:  "NUMBER  OF  DATA  TO  CHARACTERIZE  PLASMA  DISCHARGE" 

(15)  ONE  CARD. 

NPLAS  -  NUMBER  OF  SETS  OF  DATA  IN  FILE  10.0.  NPLAS  MOST  NOT 
BE  LESS  THAN  TWO  OR  GREATER  THAN  ONE  HUNDRED.  EXCEPT 
IN  THE  SPECIAL  CASE  WHEN  PBSK  (FILE  4.0)  >  0.  AND  THE 
VALUE  NPLAS  -  0  IS  THEN  ALLOWED. 

FILE  9.0:  "TOTAL  PLASMA  ENERGY  AND  MASS"  (2F10.0)  ONE  CARD. 

TOTPE  -  TOTAL  PLASMA  ENERGY  (J). 

TOTPM  -  TOTAL  PLASMA  MASS  (G) . 
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FILE  10.0:  -PIASMA  DISCHARGE  CHARACTERISTICS* 

(3F10.0)  NPLAS  CARDS. 

N.B.  THIS  FILE  IS  REQUIRED  IF  AMD  ONLY  IF  NPLAS  >  I. 

PLAST(l)  -  FIRST  VALUE  OF  TINE  (MSEC). 

PLASE(I)  -  FRACTION  OF  TOTAL  PLASMA  ENERGY  DELIVERED  AT  TINE 
PLAST(l)  (-). 

PLASM(l)  -  FRACTION  OF  TOTAL  PLASMA  MASS  DELIVERED  AT  TIME 
PLAST(l)  (-). 

PLAST(2)  -  SECOND  VALUE  OF  TIME  (MSEC) . 


PLASM(NPLASM) 


NOTES:  AT  TIMES  OUTSIDE  THE  TABLE  RANGE  THE  FIRST  OR  LAST 

ENTRIES  ARE  USED.  UMEAR  INTERPOLATION  IS  USED 
INSIDE  THE  TABLE  RANGE. 

PLASE  AND  PLASM  ARE  EXPECTED  TO  BE  NONOTONIC 
FUNCTIONS  OF  TIME  AMD  TO  BE  GREATER  THAN  OR  EQUAL 
TO  ZERO  AMD  LESS  THAN  OR  EQUAL  TO  ONE. 

IF  PBSK  (FILE  4.0)  IS  MOT  ZERO.  THE  VALUES  OF 
PLASE  AND  PLASM  ARE  NOT  USED. 

THE  TABLES  SHOULD  BE  SUFFICIENTLY  COMPLETE  TO  SPECIFY 
THE  CONDITION  AT  WHICH  THE  EFFECTIVE  ENERGY  OF  THE 
MIXTURE  IS  ZERO.  OTHERWISE  CALCULATIONS  IN  WHICH  THE 
PLASMA  DISCHARGE  YIELDS  ENERGY  INPUT  BELOW  THE  TABLE 
RANGE  WILL  EXHIBIT  ERRONEOUSLY  HIGH  PRESSURES  DUE 
TO  THE  DEFAULT  PROCEDURES. 

FILE  11.0:  "NUMBER  OF  DATA  TO  CHARACTERIZE  RATE  OF  DECOMPOSITION 
OF  WORKING  FLUID*  (15)  ONE  CARD. 

NALPH  •  NUMBER  OF  SETS  OF  DATA  IN  FILE  12.0.  NALPH  MUST  BE 

AT  LEAST  TWO  AND  NOT  GREATER  THAN  ONE  HUNDRED  EXCEPT 
IN  THE  SPECIAL  CASES  WHEN  RATE  (FILE  4.0)  >  0.  OR 
PBSK  >  0  (FILE  4.0).  OR  NBAL  (FILE  15.0)  >  0.  AMD  THE 
VALUE  NALPH  -  0  IS  THEN  ALLOWED. 


FILE  12.0:  "RATE  OF  DECOMPOSITION  OF  WORKING  FLUID* 

(2F10.0)  NALPH  CARDS. 

N.B.  THIS  FILE  IS  REQUIRED  IF  AND  ONLY  IF  NALPH  >  1. 

ALPT(l)  -  FIRST  VALUE  OF  TIME  (MSEC). 

ALPH(l)  -  FRACTION  OF  WORKING  FLUID  DEC(»(POSED  AT  TIME  ALPT(l) 

(-). 

ALPT(2)  -  SECOND  VALUE  OF  TIME.  STARTS  A  NEW  CARO. 


ALPH(NALPH) 

NOTES:  IF  RATE  (FILE  4.0)  IS  NOT  ZERO,  OR  IF  NBAL  (FILE  15.0) 

IS  NOT  ZERO,  THE  DATA  IN  THIS  FILE  ARE  NOT  USED. 

THE  VALUES  OF  ALPH  ARE  EXPECTED  TO  BE  A  MONOTONIC 
FUNCTION  OF  TIME  AND  TO  BE  GREATER  THAN  OR  EQUAL  TO 
ZERO  AND  LESS  THAN  OR  EQUAL  TO  ONE. 
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FILE  13.0:  "NUMBER  OF  DATA  TO  CHARACTERIZE  COMPOSITION  OF 

DEC(»(FOSED  WORKING  FLUID  AS  A  FUNCTION  OF  PLASMA 
ENERGY"  (15)  ONE  CARD. 

NVDAT  -  NUMBER  OF  SETS  OF  DATA  IN  FILE  14.0.  IF  NVDAT 
IS  ENTERED  AS  ZERO,  THE  COMPOSITION  OF  THE 
PROPULSION  GASES  IS  ASSUMED  TO  BE  CONSTANT  AND 
TO  BE  CHARACTERIZED  BY  THE  DATA  OF  FILE  4.0. 
NVDAT  MAY  NOT  BE  GREATER  THAN  100. 


FILE  14.0; 


N.B.  THIS 
VEIN(l) 

VEFF(l) 

VGAM(l) 

VB(1) 

VMOL(l) 

VECHB(l) 


VEIN(2) 


"PROPERTIES  OF  FINAL  DECOMPOSITION  OF  WORKING  FLUID 
AS  A  FUNCTION  OF  TOTAL  INJECTED  PLASMA  ENERGY" 
(6F10.0)  NVDAT  CARDS. 

FILE  IS  REQUIRED  IF  AND  ONLY  IF  NVDAT  IS  NOT  ZERO. 

-  FIRST  VALUE  OF  PLASMA  ENERGY  PER  UNIT  MASS 
OF  WORKING  FLUID  (J/G) . 

-  CORRESPONDING  VALUE  OF  EFFECTIVE  INTERNAL  ENERGY 
(J/G) . 

-  CORRESPONDING  VALUE  OF  RATIO  OF  SPECIFIC  HEATS  (-) 

-  CORRESPONDING  VALUE  OF  COVOLUME  (CC/G) . 

-  CORRESPONDING  VALUE  OF  MOLECULAR  WEIGHT  (G/GMOL) . 

-  CORRESPONDING  VALUE  OF  DIFFERENCE  IN  HEATS  OF 
FORMATION  OF  INTERMEDIATE  AMD  FINAL  PRODUCTS  (J/G) 
THIS  QUANTITY  IS  TAKEN  TO  BE  POSITIVE  WHEN  THE 
REACTION  FROM  INTERMEDIATE  TO  FINAL  PRODUCTS  IS 
EXOTHERMIC. 

•  SECOND  VALUE  OF  PLASMA  ENERGY  PER  UNIT  MASS 
OF  WORKING  FLUID  (J/G).  STARTS  A  NEW  CARD. 


VECHB (NVDAT) 


NOTES:  IF  THE  ENERGY  DENSITY  OF  THE  MIXTURE  OF  PLASMA 

AND  DECCMfPOSED  WORKING  FLUID  LIES  OUTSIDE  THE 
RANGE  OF  THE  TABLE  THE  FIRST  OR  LAST  VALUES  ARE 
USED  TO  DETERMINE  THE  EFFECTIVE  PROPERTIES  OF 
THE  MIXTURE. 

THE  ONLY  EXCEPTION  TO  THIS  RULE  OCCURS  FOR  THE 
EFFECTIVE  INTERNAL  ENERGY  AT  VALUES  OF  ENERGY 
DENSITY  LARGER  THAN  VEIN (NVDAT) .  THE  EFFECTIVE 
INTERNAL  ENERGY  IS  DETERMINED  IN  THIS  CASE  BY 
LINEAR  EXTRAPOLATION  OF  THE  LAST  TWO  VALUES  IN 
THE  TABLE. 

INSIDE  THE  TABLE  RANGE  LINEAR  INTERPOLATION  IS 
USED  TO  DETERMINE  VALUES  OF  ALL  PROPERTIES. 
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FILE  15:  "NUMBER  OF  SETS  OF  DATA  TO  DEFINE  INTERIOR  BALLISTIC 
BEHAVIOR  OF  GUN"  (IS)  ONE  CARD. 

NBAL  -  NUMBER  OF  SETS  OF  DATA  IN  FILE  16.  NBAL  MAY  BE 

EITHER  ZERO  OR  ANY  NUMBER  GREATER  THAN  OR  EQUAL 
TO  TWO  AND  LESS  THAN  OR  EQUAL  TO  ONE  THOUSAND. 

IF  NBAL  IS  NOT  EQUAL  TO  ZERO.  AN 
INVERSE  ANALYSIS  IS  PERFORMED  WITH  THE  RATE  OF 
DECOMPOSITION  OF  THE  WORKING  FLUID  BEING 
DETERMINED  FROM  THE  OBSERVED  HISTORY  OF  BREECH 
PRESSURE  AND  (IF  APPLICABLE)  THE  HISTORY  OF 
BASE  PRESSURE.  TOGETHER  WITH  A  CHARACTERIZATION 
OF  THE  RATE  OF  PIASMA  DISCHARGE  AND  THE  THERMO- 
CHEMICAL  PROPERTIES  OF  THE  MIXTURE. 

NSMTH  -  IF  NSMTH  IS  ZERO.  THE  BALLISTIC  DATA  ARE  NOT  SMOOTHED. 

IF  NSMTH  IS  ONE.  THE  DATA  ARE  SMOOTHED  USING  A  HIGH 
FREQUENCY  NUMERICAL  FILTER. 

NDTR  -  NUMBER  OF  DATA  TO  BE  USED  IN  REGRESSION  ANALYSIS  OF 

RATE  OF  DECOMPOSITION  OF  WORKING  FLUID.  SUGGESTED 
VALUE  IS  10.  NDTR  MAY  TAKE  ANY  VALUE  BETWEEN  ZERO  AND 
TWENTY.  IF  NDTR  IS  LESS  THAN  TWO,  A  REGRESSION  ANAL¬ 
YSIS  IS  NOT  USED  TO  ESTIMATE  THE  TREND  OF  THE  RATE  OF 
DECOMPOSITION. 

FILE  16:  "OBSERVED  BALLISTIC  RESPONSE  OF  GUN"  (3F10.0)  NBAL  CARDS. 

N.B.  THIS  FILE  IS  REQUIRED  IF  AND  ONLY  IF  NBAL  >  1. 

BALT(l)  -  FIRST  VALUE  OF  TIME  (MSEC). 

BALP(l)  -  CORRESPONDING  VALUE  OF  BREECH  PRESSURE  (MPA). 

BALPB(l)  -  CORRESPONDING  VALUE  OF  BASE  PRESSURE  (MPA). 

BALT(2)  -  SECOND  VALUE  OF  TIME.  STARTS  A  NEW  CARD. 


BALPB(NBAL) 

NOTES:  IF  VALUES  OF  BALPB  ARE  NOT  AVAILABLE  THE  ANALYSIS 

WILL  USE  THE  VALUE  OF  BET  GIVEN  IN  FILE  7.0  TO 
CHARACTERIZE  THE  MOTION  OF  THE  UNDECOMPOSED 
WORKING  FLUID.  IF  VALUES  OF  BALPB  ARE  GIVEN.  THE 
ANALYSIS  WILL  DETERMINE  AN  EFFECTIVE  VALUE  OF  BET 
AT  EACH  TIME  STEP  IN  ADDITION  TO  THE  RATE  OF 
DECOMPOSITION  OF  THE  WORKING  FLUID. 


PARAMETER  (NYI-4,NYPD-4) 

CHARACTER*10  FILNAM.OUTFIL 
CHARACTER*80  TITLE 

COMMON/Cl/  PRM,PRTRV,DB,AB,VO,C,RHOPO,AK1,AK2,QV,B,GAM,BET, 

*  DTNOM , XMO , PO , G , EGO , GMOL , TEMP , PBR , ALF . EFRAC , RATE , 

*  PBSK , PI , RU , ALFOLD , TOTPE , TOTPM , PIM , PLE 

*  .GAMO.BO.GMOLO.ECHBO.EGOO 
COMMON/C2/  RESX(20) ,RESP(20) .NRES.NSHK 

COMMON/C3/  PIAST(IOO) .PLASE(IOO) .PLASM(IOO) .NPLAS.MORE 
C0MM0N/C4/  ALPT(IOO) ,ALPH(100) .NALPH 
COMMON/C 5/  TITLE 
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C0MM0N/C6/  TIME,YI(Nyi).YDP(NYPD),YDT(NYI) 

C0MM0N/C7/  VEIN(IOO) ,VEFF(100) .VGAM(IOO) .VB(IOO) .VMOL(IOO) , 

*  VECHB(IOO) .NVDAT 

CONNON/C8/  ARCS, ARXS. AREB. AROB, ECHB. BI ,GANI ,GNOLI . CVI ,CPI 
COMMON/C9/  BALT(IOOO) .BALP(IOOO) .BALPB(IOOO) .NBAL 
COMMON/CIO/  ZA(20).RA(20),NSTA.HSPLT 
COMMON/Cll/  DXPLE.DXPIM.OXALF.RHO.VOL 
COMMON/C12/  DAT(20) .DATA(20) .DATEX.AEX.BEX.NDTR.NSTEP 
COMMON/C13/  PTOL.PTOLC.TTOL,ATOLl,ATOL2 
C 

EQUIVALENCE  (YI(1)  ,XP)  .  (YI(2)  .VP)  ,  (YI(3)  ,WF)  .  v'YI(4)  ,YFR)  , 

*  (YDP(l) ,PC) . (YDP(2) ,PB) , (YDP(3) ,RHOP) , (YDP(4) ,RES) . 

*  (YDT(2).VPRDOT) 

C 

C  NAME  FILES  FOR  INPUT  AND  OUTPUT  ON  MICROCOMPUTER 
C 

C*MFST  WRITE(*,1) 

C*MFST  1  FORMATC  ENTER  INPUT  DATA  FILE  NAME:  '\) 

C*MFST  READ(*,2)  FILNAM 

C*MFST  2  FORMAT (AlO) 

C*MFST  OPEN(5,FILE-FILNAM) 

C*MFST  WRITE(*,3) 

C*MFST  3  FORMATC  ENTER  OUTPUT  FILE  NAME:  '\) 

C*MFST  READ(*,2)  OUTFIL 

C*MFST  OPEN(6,FILE-OUTFIL) 

C 

C  READ  INPUT  DATA 
C 

READ(5.1004)  TITLE 

READ(5,1000)  DTNON,PTOL,PTOLC,TTOL,ATOLl,ATOL2 
IF(PTOL.LT.l.E-lO)  PT0I/-1.E-4 
IF(PTOLC.LT.l.E-lO)  PTOLC-l.E-4 
IF(TTOL.LT.l.E-lO)  TTOL^l.E-4 
IF(ATOLl.LT.l.E-lO)  ATOLl-l.E-4 
IF(AT0L2.LT.1.E-10)  ATOL2-1.E-2 
READ (5. 1001)  NPRT 

READ( 5 , 1000)  PRM , PRTRV , DB , VO , TEMPO , PO , PBSK , RATE 
READ(S,1001)  NSTA.NSPLT 

IF(NSTA.GE.2)  READ(5,1002)  (ZA(I) ,RA(I) , I-l ,NSTA) 

READ( 5,1001)  NRES.NSHK 

IF(NRES.GT.O)  READ(5,1002)  (RESX(I) ,RESP(I) , I-1,NRES) 
R£AO(5,1000)  C.RHOPO.AK1,AK2,QV,BET 
READ(5.1000)  B.GAM.GMOL 

REA0(S,1006)  ARCB.ARXB.AREB.AROB.ECHB.BI.GAMI.GMOLI 
READ (5, 1001)  NPLAS 
READ(5,1000)  TOTPE.TOTPM 
IF(NPLAS.GE.2) 

*  READ(5,1003)  (PLAST(I) ,PLASE(I) ,PLASM(I) ,1-1, NPLAS) 

READ (5, 1001)  NALPH 

IF(NALPH.GE.2)  READ(5,1002)  (ALPT(I) ,ALPH(I) , I-l, NALPH) 
READ(5,1001)  NVDAT 

IF(NVDAT.GE.2)  READ(5 . 1005) (VEIN(I) ,VEFF(I) ,VGAM(I) , VB(I) , 

*  VM0L(I),VECHB(I),I-1,NVDAT) 
READ(5,1001)  NBAL,NSMTH,NDTR 
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IF(NRAL.GE.2)  READ(S.I003)  (BALT(I) ,BALP(I) .BALPB(I) . I-l.NBAL) 


WRITE  INPUT  DATA 
WRITE(6,2000)  TITLE 

WRITE(6 . 2010)  DTNOM. PTOL, PTOLC.TTOL, ATOLl . AT0L2 
WRITE(6.2020)  MFRT 

WRITE(6 . 2030)  PRM, PRTRV . DB .VO .TEMPO . PO . PBSK . RATE 
WRITE(6.2032)  NSTA.MSPLT 

IF(NSTA.GE.2)  WRITE(6.2034)  (ZA(I) .RA(I) . I-l.NSTA) 

WRITE(6.2040)  NRES.NSHK 

IF(NRES.GT.O)  URITE(6.20S0)  (RESX(I) .RESP(I) . I-l ,NRES) 
WRITE(6.2060)  C.RHOPO.AK1.AK2.QV.BET 
WRITE(6.2062)  B.GAM.GMOL 

WRITE(6.2064)  ARCB.ARXB.AREB.AROB.ECHB.BI.GAMI .GMOLI 
WRITE(6.2070)  NPLAS 
WRITE(6.2080)  TOTPE.TOTPM 
IF(NPLAS.GE.2) 

*  WRITE(6.2090)  <PLAST(I).PLASE(I).PLASM(I) ,1-1, NPLAS) 
WRITE(6.2100)  NALPH 

IF(NALPH.GE.2)  WRITE(6 , 2110)  (ALPT(I) ,ALPH(I) . I-l .NALPH) 
WRITE(6.2120)  NVDAT 

IF(NVDAT.GE.2)  WRITE(6 . 2130)  (VEIN(I) , VEFF(I) , VGAM(l) . VB(I) . 

*  VMOL(I),VECHB(I),I-l.NVDAT) 
WRITE(6.2140)  NBAL.NSMTH.NDTR 

IF(NBAL,GE.2)  WRITE(6,2150)  (BALT(l) ,BALP(I) ,BALPB(I) , I-1,NBAL) 

INITIALIZE 

PI-3.1415927E0 

G-l.E+7 

RU-8.31434E0 

AB-PI*DB**2/4. EO 

DTNOM-DTNOM*l.E-3 

GAMO-GAM 

BO-B 

GNOLO-GHOL 

ECHBO-ECHB 

RESCALE  PLASMA  DATA 

IF(NPLAS.EQ.O)  GO  TO  102 
DO  100  I-l, NPLAS 
PLAST(I)-PLAST(I)*l.E-3 
PLASE ( I ) -PLASE ( I ) *TOTPE 
PLASM ( I ) -PLASM ( I ) *TOTPM 
100  CONTINUE 

102  IF(NALPH.EQ.O)  GO  TO  112 
DO  110  1-1, NALPH 
ALPT(I)-ALPT(I)*1 . E- 3 
110  CONTINUE 

112  IF(NBAL.EQ.O)  GO  TO  160 
DO  120  I-l.NBAL 
BALT( I )-BALT( I )*1 . E- 3 
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120  COHTINUE 

IF(NSimi.EQ.l)  CALL  SMOOTH 
160  XF-O.EO 
VP-O.EO 
WF-O.EO 
YFR-O.EO 
ALFOLD-O.EO 
TIME-O.EO 
FMAX-O.EO 

IF(NV1>AT.GT.0)  ECHB-VECHB(l) 

TEMP-TEMPO 

PC-PO 

CALL  EOSLIQ(PC.RHOP.DRHOP,EP,DEP) 

VFO-VO-C/RHOP 

RHOO-1 . EO/(B+RU*TEMPO/GMOL/PO) 

XMO>4(HOO«VFO 
CV-RU/GM0V(GAM-  1 .  EO) 

CP-CV*GAM 

EG00-CV*TEMP0 

EGO-EGOO^XMO-fEP 

IF(ARCB.GT.l.E-lO)  GO  TO  170 

CVI-CV 

CPI-CP 

GMOLI-GMOL 

BI-B 

GAMI-GAM 
GO  TO  180 

170  CVI-RU/GMOLI/<GAMI - 1 . EO) 

CPI-CVI*GAMI 

180  CALL  GRAD(PC,PBR,PB.XP.VP,0.E0,TOTKE,l,l,0,0) 

NSTEP-0 

MORE-1 

IF(NPRT.LT.l)  NPRT-1 
NSTOP-0 

TEST  FOR  PRINT 

200  IF(MOD(NSTEP,NPRT).NE.O)  GO  TO  210 
CALL  DEPEND 

IP(MOD(NSTEP/NPRT,50).EQ.O)  HRITE(6,3000)  TITLE 
205  TIMMS-TIME*! .  E4-3 
VPM-VP*l.E-2 

WRITE(6,3010)  TIMMS, PBR.PB, TEMP, YFR.XP.VPM.ALF, BET, EFRAC 
IF(NSTOP.BQ,0)  GO  TO  210 
HRITE(6 , 3020)  PMAX 
STOP 

TEST  FOR  MAXIMUM  PRESSURE 
210  IF(PBR.GT.PMAX)  PMAX-PBR 

SET  TIME  STEP  AMD  PREPARE  FOR  STOP 


DT-DTNOM 
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IF(XP-»-VP*DT+VPRDOT*DT*DT*0.5E0.LT.PRTRV)  GO  TO  240 
NSTOP^l 

IF(  ABS(VPRDOT).LT.l.E-lO)  GO  TO  230 

DT-(  SQRT(VP*VP+2.EO*VPRDOT*(PRTRV-XP))-VP)/VPRDOT 

GO  TO  240 

230  IF(  ABS(VP).LT.l.E-lO)  GO  TO  240 
DT-(PRTRV-XP)/VP 

UPDATE  TIME  AND  INDEPENDENT  STATE  VARIABLES 

240  CALL  RKUT(DT) 

ALF0L0>^ALF 
NSTEP^STEP-M 
IF(NDTR.LT.2)  GO  TO  260 
DO  250  I-2.NDTR 
DAT(I-1)-DAT(I) 

DATA(I-1)-DATA(I) 

250  CONTINUE 

DAT (NDTR) -TIME 
DATA(MDTR)-ALF 
260  IF(NSTOP.EQ.O)  GO  TO  200 
GO  TO  205 

FORMAT  STATEMENTS 

1000  F0RMAT(8F10.0) 

1001  FORMAT(16I5) 

1002  F0RMAT(2F10.0) 

1003  F0RMAT(3F10.0) 

1004  F0RMAT(A80) 

1005  F0RMAT(6F10.0) 

1006  FORMAT (ElO. 4, 7F10.0) 

2000  FORMATC'  '  ,10X. 'LUMPED  PARAMETER  SIMULATION  OF  ELECTROTHERMAL  GUN 
*'//10X,A80//) 

2010  FORMATC 

★'  INTEGRATION  TIME  STEP(MS)  ',F10.3/ 

*'  PRESSURE  CONVERGENCE  TOLERAMCE(-)  ' .G12.4/ 

*'  CONVERGENCE  TOLERANCE  FOR  CONSTANT  PRESSURE'/ 

*'  SOLUTION(-)  ',G12.4/ 

*'  TEMPERATURE  CONVERGENCE  TOLERANCE(-)  ',G12.4/ 

*'  FIRST  CONVERGENCE  TOLERANCE  FOR  FRACTION'/ 

*'  OF  WORKING  FLUID  DECOMPOSED( - )  ',G12.4/ 

*'  SECOND  CONVERGENCE  TOLERANCE  FOR  FRACTION'/ 

*'  OF  WORKING  FLUID  DECOMFOSED(-)  ',G12.4//) 

2020  FORMATC 

*•  NUMBER  OF  STEPS  FOR  SOLUTION  PRINTOUT  '.15) 

2030  FORMATC 


*» 

PROJECTILE  MASSCG) 

' ,F10.2/ 

** 

PROJECTILE  TRAVELCCM) 

' ,F10.2/ 

** 

BORE  DIAMETERCCM) 

' ,F10.2/ 

★» 

VOLUME  OF  MIXING  CHAMBERCCC) 

' .FlO.l/ 

★  » 

INITIAL  TEMPERATURE  OF  GASCK) 

' .FlO.l/ 

★  » 

INITIAL  PRESSURE  CMPA) 

' .FlO.l/ 

*» 

CONSTANT  BREECH  PRESSURE CMPA) 

' .FlO.l/ 
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*'  RATE  OF  DECOMPOSITION  OVER  RATE  OF  ENERGY'/ 

*'  ADDITION(-)  '.F10.3//) 

2032  FORMAT( 

NUMBER  OF  DATA  TO  DESCRIBE  CHAMBER  GEOMETRY  ' ,15/ 

(0  IMPLIES  LAGRANGE  GRADIENT;  >I  IMPLIES  CHAMBRAGE  GRADIENT)'/ 
*'  NUMBER  OF  SUBDIVISIONS  FOR  CHAMBER  INTEGRALS  ',15//) 

2034  FORMAT( '  ' , 5X , ' CHAMBER  GEOMETRY' // 

*'  ',5X,'DISTANCE(CM)',5X,'RADIUS(CM)'//('  '.2F15.3)) 

2040  FORMAT( 

*•  NUMBER  OF  BORE  RESISTANCE  DATA  ' ,15/ 

*'  AIR  SHOCK  RESISTANCE  (O-NO; 1-YES)  ',15//) 

2050  FORMAT('  ' , 5X, 'TRAVEL(CM) ' , 5X. 'RESISTANCE(MPA) '// 

*('  ' .F15.2,5X,F15.1)) 

2060  FORMAT ('  '// 

*'  ' ,5X, 'PROPERTIES  OF  UNDECOMPOSED  WORKING  FLUID'// 


*'  MASS  OF  WORKING  FLUID(G)  '.FlO.l/ 
*'  DENSITY  OF  FLUID  AT  1  ATM.(G/CC)  ',F10.4/ 
*'  BULK  MODULUS  AT  1  ATM. (MPA)  FlO.l/ 
*'  DERIVATIVE  OF  MODULUS  W.R.T.  PRESSURE(-)  ',F10.3/ 
*'  HEAT  OF  VAPORIZATION(J/G)  ', FlO.l/ 


* '  LAGRANGE  COUPLING  COEFFICIENT( - >  ' , FIO . 2/ 

*'  (O-LIQUID  AT  REST; 1-LIQUID  HAS  LAGRANGE  VELOCITY)'//) 

2062  FORMATC  ', 5X, ' PROPERTIES  OF  FINAL  PRODUCTS  OF  DECOMPOSITION'// 

*'  COVOLUME  (CC/G)  ',F10.4/ 

*'  RATIO  OF  SPECIFIC  HEATS  (-)  ' ,F10.4/ 

*'  MOLECULAR  WEIGHT  (G/GMOL)  ' ,F10.4//) 

2064  FORMAT('  ' ,5X, 'PROPERTIES  OF  INTERMEDIATE  PRODUCTS  OF  DECOMPOSITIO 
*N'// 


*•  PRE-EXPONENT  IN  ARRHENIUS  REACTION  RATE  LAW'/ 


*'  (UNITS  YIELD  G/CC-SEC)  ',G15.5/ 

*'  TEMPERATURE  EXPONENT  IN  REACTION  RATE  LAW  (-)  ',G15.5/ 

*'  ACTIVATION  ENERGY  (J/GMOL)  ',G15.5/ 

*'  REACTION  ORDER  (-)  ',G15.5/ 

*'  CHEMICAL  ENERGY  OF  REACTION  (J/G)  ',G15.5/ 

*'  COVOLUME  (CC/G)  ',F10.4/ 

*'  RATIO  OF  SPECIFIC  HEATS  (-)  * ,F10.4/ 

*'  MOLECULAR  WEIGHT  (G/GMOL)  '.FIO. 4//) 

2070  FORMAT( 

*•  NUMBER  OF  DATA  TO  DESCRIBE  PLASMA  DISCHARGE  ',15) 

2080  FORMAT ( 

*'  TOTAL  PLASMA  ENERGY(J)  ',G15.6/ 

*'  TOTAL  PLASMA  MASS(G)  ',F10.3//) 


2090  FORMATC  ', 5X, 'TIME (MS) ', 5X, ' FRACTION  OF  ENERGY( -)', 5X. ' FRACTION  0 
*F  MASS(-)'//('  ' ,F15.3,5X,2F15.6)) 

2100  FORMAT ('  '// 

*'  NUMBER  OF  DATA  TO  DESCRIBE  DECOMPOSITION'/ 

*'  OF  WORKING  FLUID 

2110  FORMAT('  ', 5X, ' TIME (MS) ', 5X, ' FRACTION  OF  FLUID 
*('  ' ,5X,F15.3.5X,F15.6)) 

2120  FORMAT('  '// 

*'  NUMBER  OF  DATA  TO  DETERMINE  THERMOCHEMISTRY 


’.15//) 

DECOMPOSED(-)'// 

'.15//) 
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2130  FORHAT('  ' . 5X, ' ENERGy(J/G) ' , 5X. ' EFF.ENGY. (J/G) ' , 8X. 'GAMMA( • ) ' . 
★5X. 'COVOLUME(CC/G) ' ,12X. 'MOL.  WGT. ' 8X, 'DELTA  HOF(J/G) '// 

*('  ' ,2F15.1,5X.2F15.4.10X.F15.4,F15.1)) 

2140  FORMAT('  '// 

★'  NUMBER  OF  DATA  TO  SPECIFY  OBSERVED  BALLISTICS  ',15/ 

*'  DATA  SMOOTHING  REQUIRED  (0-NO; 1-YES)  '.15/ 

*'  NUMBER  OF  DATA  FOR  REGRESSION  ANALYSIS  OF  RATE'/ 

*'  OF  DECOMPOSITION  ',15//) 

2150  FORMAT( '  ' , 5X, 'TIME(MSEC) ' , 5X, ' BREECH  PRESS .  (MPA) ' , 5X, 

★'BASE  PRESS.  (MPA)'//('  ',3G20.6)) 

3000  FORMAT( ' 1 ' . lOX , A80/ 

★'  ' .5X.'TIME(MS)' .5X,'P-BR(MPA)' .5X.'P-BASE(MPA)' . 5X, 'TEMP(K) ' , 
*4X,'YFR(-)' .3X.'PR.  DISP(CM) ' ,2X, 'PR.  VEL(M/S) ' ,2X. 'ALPHA(') ' ,4X, 
*'BETA(-)' .4X,'E-FRAC( ')'//) 

3010  FORMAT('  ' .FIO . 3 . 3F15 . 2 . FIO. 3 ,F12. 2 .F12, 1. F12 . 3 , Fll . 3 . F13 . 3) 

3020  FORMAT('  '//'  MAXIMUM  PRESSURE(MPA) :  '.FlO.l) 

END 
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SUBROUTINE  RKUT(DT) 


FOURTH  ORDER  RUNGE-KUTTA  ALGORITHM 
PARAMETER  (MYI-4 , NYFD-4) 

COMMON/C6/  TIME,YI(NYI),YDP(NYPD).YDT(NYI) 

EQUIVALENCE  (YI(1)  .XP)  ,  (YI(2)  .VP)  .  (YI(3)  .WF)  ,  (YI(4)  .YFR)  . 

*  (YDP(l) .PC) . (YDP(2) .PB) . (YDP(3) ,RHOP) . (YDP(4) ,RES) . 

*  (YDT(2) .VPRDOT) 

DIMENSION  YB(NYI).CW(4),Q(NYI,4) 

C 

DATA  CW/0.E0.0.5E0.0.5EO.1.0E0/ 

C 

TB-TIME 

DO  10  J-l.NYI 

YB(J)-YI(J) 

10  CC^INUE 
DO  100  1-1.4 
TIME-TB+CW(I)*DT 
DO  20  J-l.NYI 
YI(J)-YB(J) 

IF(I.BQ.l)  GO  TO  20 

YI ( J ) -YB ( J ) +CW( I ) *DT*Q( J , I - 1) 

IF(J.NE.4)  GO  TO  20 
IF(YI(J).GT.1.E0)  YI(J)-1.E0 
IF(YI(J).GE,0.E0)  GO  TO  20 
Q(J.I-l)— YB(J)/CW(I)/DT 
YI(J)-O.EO 
20  CONTINUE 
CALL  DEPEND 
CALL  DERIV 
DO  40  J-l.NYI 
Q(J,I)-YDT(J) 

40  CONTINUE 
100  CONTINUE 
TIME-TB+DT 
DO  200  J-l.NYI 

YI (J)-YB(J)+DT/6 . E0*(Q(J , l)+2 . E0*Q(J , 2)+2 . E0*Q( j , 3)4Q(J ,4) ) 

200  CONTINUE 

IF(YFR.LT.O.EO)  YFR-O.EO 
IF(YFR.GT.l.EO)  YFR-l.EO 
RETURN 
END 
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SUBROUTINE  DERIV 


TIME  DERIVATIVES  OF  INDEPENDENT  STATE  VARIABLES 

PARAMETER  (NYI-4,NYPD-4) 

C 

COMMON/Cl/  PRM, PRTRV, DB . AB , VO . C .RHOPO .AKl , AK2 .QV . B .GAM, BET . 

*  DTNOM.XMO .  PO  ,G .  EGO  .OfOL,  TEMP.  PBR.ALF,  EFRAC  .RATE , 

*  PBSK,PI,RU.ALPOLD.TOTPE,TOTPM,P1M.PLE 

*  .GAMO.BO.GMOLO.ECHBO.EGOO 
COMMON/C6/  TIME , YI (NYI ) , YDP(NYPD) . YDT(NYI ) 

COMMON/C8/  ARCB . ARXB . AREB , AROB , ECHB . BI . GAMI , GMOLI . CVI . CPI 
COMMOM/C9/  BALT(IOOO) .BALP(IOOO) .BALPB(IOOO) .NBAL 
COMMON/Cll/  DXPLE.DXPLM.DXALF.RHO.VOL 
C 

EQUIVALENCE  (YI(1) .XP) . (YI(2) .VP) . (YI(3) ,WF) . (YI(4) .YFR) , 

*  (YDP(l) ,PC) , (YDP(2) .PB) , (YDP(3) .RHOP) . (YDP(4) .RES) , 

*  <YDT(2).VPRDOT) 

C 

YDT(1)-VP 

YDT ( 2 ) -G*AB/PRM* ( PB - RES ) 

IF(YDT(2).LT.O.EO)  YDT(2)-O.EO 

YDT(3)-AB*RES*VP 

YDT(4)-0.E0 

IF(ARCB.LT.l.E-lO)  RETURN 

BIT-AREB/RU/TEMP 

EXPON-O.EO 

IF(BIT.LT.50.E0)  EXPON-  EXP(-BIT) 

YDT(4)-(  (1 .  EO-YFR)*C*DXALF-YFR*DXPLM)/RHO/VOL 

*  -ARCB/RHO*TENP**ARXB*(YF11*RHO)**aROB*EXPON 
RETURN 

END 
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SUBROUTINE  DEPEND 

CALCUIATION  OF  DEPENDENT  STATE  VARIABLES 
PARAMETER  (NYI-4,NYPD-4) 

COMNON/Cl/  PRM.PRTRV,DB,AB.VO.C.RHOPO,AK1,AK2,QV.B.GAM.BET. 

*  DTN(X<.XMO,PO,G.  EGO, QfOL.TEMP.PBR.ALF.EFRAC. RATE. 

*  PBSK.PI.RU.ALPOLD.TOTPE.TOTPM.PLM.PLE 

*  .GAMO.BO.GMOLO.ECHBO.EGOO 
C(M1MON/C2/  RESX(20).RESP(20),NR£S.NSHK 

C(»<MON/C3/  PLAST(IOO) .PLASE(IOO) .PLASM(IOO) .NPLAS.MORE 
C0MM0N/C4/  ALPT(IOO) .ALPH(IOO) .NALPH 
COMMON/C6/  TIME,YI(NYI).YDP(NYPD),YDT(NYI) 

COMMON/C7/  VEIN(IOO) ,VEFT(100) ,VGAM(100) ,VB(100) .VMOL(IOO) , 

*  VECHB(100),NVDAT 

COMMON/G8/  ARCB . ARXB . AREB . AROB , ECHB , BI . GAMI , GMOLI . CVI . CPI 
COMMON/C9/  BALT< 1000) , BALP( 1000) , BALPB( 1000) , NBAL 
COMMON/Cll/  DXPLE.DXPLM.DXALF.RHO.VOL 
C0MM0N/C12/  DAT(20) ,DATA(20) .DATEX.AEX.BEX.NDTR.NSTEP 
C0MM0N/C13/  PT0L,PT0LC,TT0L.AT0U,AT0L2 

EQUIVALENCE  (YI(1) ,XP) , (YI(2) .VP) , (YI(3) ,WF) , (YI(4) ,YFR) , 

*  (YDP(l) .PC) . (YDP(2) ,PB) , (YDP(3) ,RHOP) , (YDP(4) ,RES) , 

*  (YDT(2) .VPRDOT) 

DATA  AIRGAM/1 . 4E0/ , AIRPO/O . lOlEO/ , AIRTO/293 . EO/ , 

*  AIRMtf/28.96E0/ 


QW-O.EO 

RES-O.EO 

IF(NRES.EQ.O)  GO  TO  S 

CALL  INTERP(XP,RES,DXRES.RESX.RESP,O.NRES) 

5  IF(NSHK.EQ.O)  GO  TO  10 

RESISTANCE  DUE  TO  SHOCK  AHEAD  OF  PROJECTILE 

BIT-(AIRGAM- 1 . EO) /(AIRGAM+l . EO) 

BITO  SQRT(G*AIRGAM*RU*AIRTO/AIRMV) 

BITl^l.EO-BIT 

BITU-(VP+  SQRT(VP*VP+4 . EO*BITL*BITL*BITC*BITC) ) 

*  /(2.E0*BITL*BITC) 

RES-RES+AIRPO*( (1 . EO+BIT)*BITD*BITD-BIT) 

10  IF(PBSK.GT.1,E-10)  GO  TO  20 

PLASMA  ENERGY  FROM  TABLE 

CALL  INTERP(TIME,PLE,DXPLE.PLAST,FLASE,0,MPLAS) 
CALL  INTERP ( TIME , PLM , DXPLM . PLAST , PLASM . 0 , NPLAS ) 
EFRAOPLE/PLASECNPLAS  ) 

IF(NBAL.EQ.O)  GO  TO  14 

RATE  OF  DECOMPOSITION  TO  MATCH  OBSERVED  BALLISTICS 


NPASS-2 

ITER-0 

NSRCH-0 

ALF-ALFOID 

C 

C  COMPUTE  EXPECTED  VALUE  FROM  LINEAR  REGRESSION  LINE 
C 

OATEX-O.EO 

IF(MDTR.GT.I.AND.NSTEP.GE.NDTR)  CALL  REGR 

IF(DATEX.GT.l.E-lO)  ALF-DATEX 

ALFl-ALF 

CALL  INTERP( TIME. PBRSK.DXPBR. BALT. BALP.O, NEAL) 

IF(BALPB(NBAL).LT.1.E-10)  GO  TO  13 

CALL  INTERPCTIME.PBASK.BITAXL.BALT.BALPB.O.NBAL) 

C 

C  CALL  TO  GRAD  TO  SET  SLIP  COEFFICIENT 
C 

12  IF(BALPB(NBAL).LT.1.E-10)  GO  TO  13 

CALL  GRAD(PSEEK.PBRSK.PBASK,XP.VP,RES,TOTKE.2.0.0.1) 

C 

C  CALL  TO  GRAD  TO  GET  SPACEMEAN  PRESSURE  AND  KINETIC  ENERGY 
C 

13  CALL  GRAD(PSEEK.PBRSK.PB.XP.VP.RES.TOTKE,2.0.0.0) 

IF(TIME . LT . 1 . E - 6 , AND . PSEEK . LT . PO)  PSEEK-PO 

GO  TO  28 

14  IF(RATE.GT.l.E-lO)  GO  TO  16 
C 

C  RATE  OF  DECOMPOSITION  FROM  TABLE 
C 

CALL  INTERP(TIME.ALF.DXALF,ALFT.ALPH,0,NALPH) 

GO  TO  18 
C 

C  RATE  OF  DECOMPOSITION  PROPORTIONAL  TO  PLASMA  ENERGY  FLUX 
C  ACCORDING  TO  VALUE  OF  RATE. 

C 

16  ALF-EFRAC 

IF(TIME.GT. PLAST(NPLAS) )  ALF-TIME/PLAST(NPLAS) 
ALF-ALF*RATE 
IF(ALF.GT.l.EO)  ALF-l.EO 
18  NPASS-2 
C 

C  CALL  GRAD  TO  GET  KINETIC  ENERGY 
C 

CALL  GRADCPC.PBR.PB.XP.VP.RES.TOTKE.l.O.O.O) 

GO  TO  28 
C 

C  PLASMA  ENERGY  TO  PROVIDE  CONSTANT  BREECH  PRESSURE  PBSK. 

C  OF  DECOMPOSITION  PROPORTIONAL  TO  ENERGY  FLUX. 

C 

20  ITER-0 
NPASS-l 
PLELf^. 

PLEH-TOTPE 

IF(MORE.EQ.O)  PLEl^PLEH 


RATE 


92 


26  PLE-0.5E0*(PLEL+PLEH) 

EFRAC-PLE/TOTPE 

PLM-TOTPM*EFRAC 

ALF-EFRAC 

IF(MORE . EQ . 0)  ALF-TIME*ERATE*RATE 
IF(ALF.GT.l.EO)  ALF-l.EO 
C 

C  CALL  GRAD  TO  GET  KINETIC  ENERGY 
C 

CALL  GRAO(PSEEK.PBSK,PB.XP.VP.RES,TOTKE,2.0.0.0) 

28  NIT-0 

30  CALL  EOSLIQ(PC.RHOP,DRHOP.EP.DEP) 

VOL-VO+AB*XP 

YFRM-l.EO-YFR 

BITO-PLM+XMO+ALF*C 

BIT2-V0L-B*BIT0 

BIT3— (l.E0-ALF)*C 

BIT4-BIT2+BIT3/RHOP 

BIT5-VOL+BIT3/RHOP 

RHO-BITO/BIT5 

DRHO-BITO/BIT5/BIT5*BIT3/RHOP/RHOP*DRHOP 

NITE-0 

32  PLEIN-PLE+EGO-QW-WF-PRM/G/2 . E0*VP**2-T0TKE 
*  .ALF*C*YFR*ECMB-EP 

IF(YFR. LT . . 999E0)  PLEIN-PLEIN-BITO*YFR*CVI*TEMP 
IF(NVDAT.EQ.O.OR.YFR.GE. .999E0)  PLEIN-PLEIN+ALP*C*(ECHB-QV) 
PLEIN-PLEIN/BITO 

IF(YFR . LT . . 999E0)  PLEIN-PLEIN/YFRM 
IF(NVDAT.GT.1.AND.YFR.LT. .999E0)  GO  TO  36 
C 

C  FIXED  THERMOCHEMISTRY 
C 

EMIX-PLEIN 
DXEMIX-l.EO 
DXECHB-O.EO 
DXTEMP-O.EO 
DPLEIN-O.EO 
DXGAM-O.EO 
DXGMOL-0 . EO 
GO  TO  40 
C 

C  VARIABLE  THERMOCHEMISTRY 
C 

36  CALL  INTERP(PLEIN.EMIX,DXEMIX.VEIN,VEFF.1.NVDAT) 
IF(EMIX.GE.EGOO)  GO  TO  38 
EMIX-EGOO 
GAM-GAMO 
B-BO 

gmolh:molo 

ECHB-ECHBO 
GO  TO  40 

38  CALL  INTERP(PLEIN. GAM. DXGAM, VEIN, VGAM.O.NVDAT) 

CALL  INTERP(PLEIN,B,DXB,VEIN,VB,0,NVDAT) 

CALL  INTERP ( PLEIN , GMOL . DXGMOL . VEIN . VMOL . 0 . NVDAT) 


93 


O  O  f>  o  o  o 


CALL  INTERP(PLE1H , ECHB . DXECHB . VEIN . VECHB , 0 , NVDAT) 

40  CV-RU/GMOL/(GAM-1.EO) 

CP-GAM*CV 

CV-YFRM*CV+YFR*CVI 

CP-YFRM*CP+YFR*CPI 

B-YFRM*B+YFR*B1 

GMOLF-GMOL 

GMOI^l . EO/(YFRM/GMOL+YFR/GMOLI ) 

GAM-CP/CV 

IF(YFR.GT.1.E-3.AND.YFR.LT.0.999E0)  GO  TO  42 

TEMP-EMIX/CV 

GO  TO  48 

NEED  TO  ITERATE  FOR  VALUE  OF  TEMPERATURE 

42  FPE-TEMP-EMIX/CV 

IF(  ABS(FPE) .LT.TTOL*TEMP)  GO  TO  48 
IF(NITE.LT.50)  GO  TO  46 
WRITE(6,44)  FPE.TEMP 

44  FORMAT  ('  MORE  THAN  30  ITERATIONS  TO  DETERMINE  TEMPERATURE  IN  DEPEN 
*D.  TERMINATING.'/'  FPE  ,G15.6. '  TEMP  ,G15 . 6) 

STOP 

46  DXCV-YFRM*RU*(DXGAM/GMOL-  (GAM- 1 .  EO)/GMOLF/GKOLP*DXGMOL) 

DFPE-1 . EO+(DXEMIX/CV-EMIX/CV/CV*DXCV)*YFR*CVI/YFRM 
TEMP-TEMP - FPE/DFPE 
NITE-NITE+l 
GO  TO  32 

48  IF(YFR.LT, .999E0)  DXTEMP-RU/GMOL*((l.EO/RHO-B) -PC/RHO/RHO*DRHO) 
BIT— YFR 

IF(NVDAI.EQ.O)  BIT-1. EO-YFR 

DPLEIN-( - DEP-BITO*YFR*CVI*DXTEMP)/BITO 

IF(YFR.LT. .999E0)  DPLEIN-DPLEIN/YFRM/(1.E0-BIT*ALF*C*DXECHB/BIT0 

*  /YFRM) 

BIT1-(GAM- 1 . E0)*BIT0*EMIX 

BIT2-VOL-B*BITO 

BIT4-BIT2+BIT3/RHOP 

IF<ALF.LT. .99999E0)  GO  TO  50 
IF(YFR.GT.1.E-10.AHD.YFR.LT. .999EO)  GO  TO  50 
PC-BIT1/BIT2 
GO  TO  100 

PRESSURE  AND  LIQUID  DENSITY  ARE  COUPLED  NONLINEARLY 

50  FP-PC-BIT1/BIT4 

IF(  ABS(FP).LT.PTOL*PC)  GO  TO  100 
IF(NIT.LT.50)  GO  TO  80 
WRITE(6.60)  FT. PC 

60  FORMAT ('  MORE  THAN  50  ITERATIONS  TO  DETERMINE  PRESSURE  IN  DEPEND. 
★TERMINATING. '/ 

*'  FP  -',G15.6,'  PC  -',G15.6) 

STOP 

80  DFP-1 . EO-DRHOP*BIT1*BIT3/(BIT2*RHOP+BIT3)**2 

*  -DPLEIN*DXEMIX*(GAM- 1 . EO)*BITO*RHOP/(BIT2*RHOP+BIT3) 
PC-PC-FP/DFP 
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NIT-NIT+1 
GO  TO  30 

100  IF(NPASS.EQ.2)  GO  TO  160 
C 

C  CHECK  WHETHER  CONSTANT  PRESSURE  SOLUTION  HAS  BEEN  OBTAINED 
C 

IF(MORE.EQ.O)  GO  TO  220 

IF(  ABS(PC-PSEEK) .LT.PTOLC*PSEEK)  GO  TO  220 
IF( ITER. LT. 100)  GO  TO  120 
HRITE(6.110)  PC.PLE.PLEL.PLEH 

110  F0RHAT('  HORE  THAN  100  ITERATIONS  IN  SEARCH  FOR  CONSTANT  PRESSURE 
★SOLUTION.  TERMINATING.'/ 

*'  PC.PLE.PLEL.PLEH  '.4G15.6) 

STOP 

120  IF(PC.GT.PSEEK)  PLEH-PLE 

IF(PLEH.GT.TOTPE)  PLEH-TOTPE 
IF<PC.LT.PSEEK)  PLEL^PLE 
IF(PLET .LT.O.EO)  PLEL-O.EO 
IF(PLEL.GT.TOTPE)  PLEI^TOTPE 
ITER-ITER+1 

IF(PLEL.LT.0.999*TOTPE)  GO  TO  26 
MORE-O 

ERATE-1 . EO/TIME 
GO  TO  26 

160  IF(NBAL.EQ.O)  GO  TO  220 
C 

C  CHECK  WHETHER  OBSERVED  BALLISTICS  ARE  MATCHED 
C 

IF(ALP0LD.GT.0.999E0)  GO  TO  220 

FP-PC/PSEEK-l.EO 

IF(  ABS(FP).LT.ATOLl)  GO  TO  200 

IF( ITER. LT. 100)  GO  TO  170 

IF(NSRCH.EQ.l)  GO  TO  210 

IF(  ABS(FP) .GT.ATOL2)  GO  TO  164 

WRITE(6,162)  TIME.FP 

162  FORMAT<'  AT  TIME  -'.G12.5.'  ALPHA  SEARCH  HAS  FP  -  '.G12.5. 

*  '  CONTINUING.') 

GO  TO  200 

164  WRITE(6.16S)  TIME.PC.PSEEK.ALF.BET 

165  FORMAT('  MORE  THAN  100  ITERATIONS  TO  MATCH  BALLISTICS.  TERMINATING 
*.'/'  TIME  -'.G12.5.'  PC  -'.G12.5.'  PSEEK  -' .G12 . 5 . '  ALPHA-'. 

★  G12.5. '  BETA  -' .G12.5) 

STOP 

170  IF(ITER.GT.O)  GO  TO  175 
172  ALFW-ALF 
FFW-FP 

ALF-ALF+0.01 
GO  TO  180 

175  DFP-(FP-FPW)/(ALF-ALFW) 

ALFW-ALF 

FPW-FP 

IF(  ABS(DFP) .LT.l.E-5)  GO  TO  172 
ALF-ALF-FP/DFP 

IF(NSRCH.EQ.l.AND.ALF.LT.ALFl)  GO  TO  205 
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IF(ALF.LT.O.EO)  ALF-ALFV/2 . EO 
IF(ALF.CT.l.EO)  ALF-(l.E0+ALFtf)/2.E0 
IFCNSRCH.EQ.O)  GO  TO  180 
IF(  ABS(ALF-AUFV) .LT.l.E-6)  GO  TO  205 
180  ITER-lTER-i-1 
GO  TO  12 

CHECK  FOR  MONOTOMICITY  OF  ALPHA 

200  IF(NVDAT.EQ.O)  GO  TO  220 

IF(ALF.LT.ALF0LD-l.E-3)  GO  TO  204 
IFCNSRCH.EQ.O)  GO  TO  220 

ACCEPT  LOWER  ROOT  IF  CLOSER  TO  EXPECTED  VALUE  AND  DOWNTREND  HAS 
NOT  BEEN  ESTABLISHED 

IF(DATEX.LT.l.E-lO)  GO  TO  220 
IF(BEX.LT.-l.E-lO)  GO  TO  220 

IF(  ABS(ALF-DATEX) .LT.  ABS(ALF2-DATEX) )  GO  TO  220 
GO  TO  210 

204  IFCNSRCH.EQ.O)  ALF2-ALF 
NSRCH-1 

205  IFCALFl.GE. .989E0)  GO  TO  210 
ALFl-ALFl+O.OlEO 
ALF-ALFl 

ITER-0 
GO  TO  12 
210  ALF-ALF2 

CAIL  TO  GRAD  TO  GET  BREECH  AND  BASE  PRESSURES 

220  CALL  GRADCPC.PBR,PB,XP,VP.RES,TOTKE,1,0,0,0) 

VOL^BIT5 

RETURN 

END 
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SUBROUTINE  EOSLIQ( PIN . RHOP , DRHOP , EP . DEP) 

C 

C  DETERMINES  DENSITY  (RHOP)  AND  INTERNAL  ENERGY  (EP)  OF  WORKING 

C  FLUID  AS  A  FUNCTION  OF  PRESSURE  (PIN)  TOGETHER  WITH  DERIVATIVES 

C  DRHOP  AND  DEP. 

C 

COMMON/Cl/  PRM.PRTRV,DB,AB,VO,C,RHOPO,AK1,AK2,QV,B,GAM,BET. 

*  DTNOM , XMO , PO , G , EGO . GMOL ,  TEMP , PAR , ALF , EFRAC , RATE , 

*  PBSK.PI.RU.ALFOLD.TOTPE.TOTPM.PLM.PLE 

*  .GAMO.BO.GMOLO.ECHBO.EGOO 
C 

P-PIN-O.IOIEO 

IF(P.GT.O.EO)  GO  TO  10 

RHOP-RHOPO 

DRHOP-0 . EO 

EP-0 . EO 

DEP-O.EO 

GO  TO  80 

10  IF(  ABS(AK2).GT.l.E-6)  GO  TO  20 
RHOP-RHOPO* ( 1 . EO+P/AKl) 

DRHOP-RHOPO/AKl 
GO  TO  40 

20  RHOP-RHOPO* ( 1 . E0+AK2/AK1*P) ** ( 1 . E0/AK2 ) 

DRHOP-RHOPO/AK1*(1 . E0+AK2/AKl*P)**( (l.E0-AK2)/AK2) 

40  EP— P/RHOP 

DEP— 1 . EO/RHOP+P/RHOP/RHOP*DRHOP 

IF(  ABS(AK2).GT.1.E-6.AND.  ABS(AK2-1.E0) .GT.1.E*6)  GO  TO  60 
EP-EP+AK1/RHOPO*ALOG ( 1 . EO+P/AKl ) 

DEP-DEP+1 . EO/RHOP/( 1 . EO+P/AKl ) 

GO  TO  80 

60  EP-EP+ ( AK1+AK2*P ) / ( AK2 - 1 . EO ) /RHOP - AKl/ ( AK2 - 1 . EO ) /RHOPO 
DEP-DEP+(AK2- (AK1+AK2*P)/RH0P*DRH0P)/(AK2-1 . EO)/RHOP 
80  EP-(1.E0-ALF)*C*EP 
DEP-( 1 . EO - ALF) *C*DEP 
RETURN 
END 
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SUBROUTINE  INTERP(X,Y,DYDX.XA,YA,NXTRAP,N) 

C 

C  INTERPOLATE  (XA.YA)  ARRAYS  TO  GET  Y(X)  AND  DERIVATIVE  DYDX. 

C  IF  X  IS  OUTSIDE  RANGE  OF  XA,  Y  IS  SET  EQUAL  TO  FIRST  OR  LAST 

C  VALUE  OF  YA  UNLESS  NXTRAP  -  1  IN  WHICH  CASE  LINEAR  EXTRAPOLATION 

C  IS  USED  FOR  X  >  XA(N) . 

C  IF  N  IS  LESS  THAN  TWO.  Y  IS  SET  EQUAL  TO  ZERO. 

C 

DIMENSION  XA(N) .YA(N) 

C 

Y-O.EO 
DYDX-O.EO 

IF(N.LT.2)  GO  TO  50 
IF(X.LT.XA(1))  GO  TO  30 
IF(X.CT.XA(N))  GO  TO  40 
DO  10  1-2, N 

IF(XA(I) .GT.X)  GO  TO  20 
10  CONTINUE 
GO  TO  40 

X  IS  IN  TABLE  RANGE.  INTERPOLATE. 

20  BIT-XA(I)-XA(I-1) 

IF(BIT.GT.l.E-lO)  DYDX-(YA(I)-YA(I-1))/BIT 
Y-YA ( I - 1 ) + ( X - XA ( I - 1 ) ) *D YDX 
GO  TO  50 

X  IS  OUTSIDE  TABLE  RANGE.  USE  FIRST  OR  LAST  VALUE  OF  YA, 

30  Y-YA(l) 

GO  TO  50 

40  IF(NXTRAP.EQ.l)  GO  TO  45 
Y-YA(N) 

GO  TO  50 
45  I-N 

GO  TO  20 
50  RETURN 
END 
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SUBROUTINE  GRAD( PC , PBRCH . PB .XP . VP , RES . TOTKE . IND , IHIT , KALC , MODE) 


C 

C  ACTION  DEPENDS  ON  VALUE  OF  INPUT  DATUM  MODE.  WHEN  MODE  -  0,  GRAD 

C  CCHfPUTES  EFFECT  OF  PRESSURE  GRADIENT.  PC,  PBRCH  AND  PB  ARE 

C  SPACEMEAN,  BREECH  AND  BASE  PRESSURES  RESPECTIVELY. 

C  XP  AND  VP  ARE  PROJECTILE  DISPLACEMENT  AND  VELOCITY  RESPECTIVELY. 

C  RES  IS  BORE  RESISTANCE  PRESSURE.  TOTKE  IS  KINETIC  ENERGY  OF 
C  MOVING  FLUID. 

C  IND  -  1.  2  OR  3  ACCORDING  AS  PC,  PBRCH  OR  PB  IS  SUPPLIED  AS  INPUT 

C  TO  ROUTINE. 

C  INIT  -  1  FOR  INITIALIZATION  CALL  (NECESSARY  IF  CHAMBER  GEC»(ETRY 

C  IS  NEWLY  ENTERED  TO  PROGRAM)  AND  0  SUBSEQUENTLY. 

C  KALC  -  0  IF  THE  CHAMBER  INTEGRALS  NEED  TO  BE  UPDATED  AT  A  STEP 

C  DURING  THE  SOLUTION  AND  1  IF  THE  VALUES  ARE  ALREADY  KNOWN. 

C 

C  INPUT  ARE  XP.  VP,  RES,  IND,  INIT,  KALC  AND  MODE  AND  ONE  OF  PC,  PBR 
C  OR  PB. 

C 

C  WHEN  MODE-1,  GRAD  COMPUTES  THE  VALUE  OF  THE  SLIP  COEFFICIENT  BET 

C  (COMMON  BLOCK  /Cl/)  TO  MATCH  THE  VALUES  OF  PBRCH  AND  PB. 

C 

C  INPUT  ARE  XP,  VP.  RES,  IND,  INIT.  KALC  AND  MODE  AND  BOTH  PBRCH 
C  AND  PB. 

C 

PARAMETER  (NYI-4,NYPD-4) 

C 

COMMON/C 1 /  PRM , PRTRV ,DB,AB,VO,C,RHOPO,AK1,AK2,QV,B,GAM,BET, 

*  DTNOM,XMO,PO,G, EGO, GMOL, TEMP, PBR, ALF.EFRAC, RATE, 

*  PBSK,PI,RU,ALF0LD,T0TPE,T0TPM,PLM,PLE 

*  ,GAMO,BO,GMOLO,ECHBO,EGOO 
COMMON/C6/  TIME,YI(NYI),YDP(NYPD),YDT(NYI) 

COMMON/CIO/  ZA(20) ,RA(20) .NSTA.NSPLT 

C 

IF(INIT.EQ.O)  GO  TO  100 

INITIALIZATION 

IF(NSTA.EQ,0)  RETURN 
IF(NSTA.LT.2.0R.NSTA.GT.20)  GO  TO  320 
IF(NSPLT.LT.l)  NSPLT-1 
VOLO-O.EO 
AJIO-O.EO 
AJ30-0.E0 
AJ40-0.E0 
Z-O.EO 

CALL  INTERP(0.E0,R,DUM,ZA,RA,0,NSTA) 

AL^PI*R*R 
VOLOL-O.EO 
AJIOL-O.EO 
AJ2L-0.E0 
DO  80  I-2,NSTA 

DZ-(ZA(I)-ZA(I-1))/  FLOAT(NSPLT) 

DO  60  J-1,NSPLT 
Z-Z+DZ 
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CALL  INTERF(Z,R,DUM,ZA,RA,0,NSTA) 

A-PI*R*R 

VOLO-VOLO+DZ*(A+AL)/2 . EO 
AJ10-AJ10+DZ*(VOL0+VOL0L)/(A+AL) 

AJ2-(VOLO/A)**2 

AJ30-^30+DZ*(A+AL)*(AJ10+AJl0L)/4 .  EO 
AJ40-AJ40-i-DZ*(A+AL)*(AJ24^AJ2L)/4  .  EO 
AL-A 

VOLOI/-VOLO 
AJlOb-AJlO 
AJ2I/-AJ2 
60  CONTINUE 
80  CONTINUE 

WRITE(6.90)  VOLO,AJ10,AJ30.AJ40 

90  FORMATC//'  GRAD  INTEGRALS  ...  VOLO: \G12.4, '  AJIO: ' ,G12.4, '  AJ30;' 
*,G12.4,*  AJ40:' .G12.4//) 

RETURN 

C(»1PUTE  PRESSURE  GRADIENT 

100  IF(MODE.EQ.O)  GO  TO  102 
IF(NSTA.BQ.O)  GO  TO  140 
GO  TO  200 

102  IF(IND.LT.1.0R.IND.GT.3)  GO  TO  300 

XK-(C*(ALF+BET*(1 . EO-ALF) )+PLM+XMO)/PRM 
IF(NSTA.GT.O)  GO  TO  200 

LAGRANGE  GRADIENT  ANALYSIS 

XK2«1.E0-^XK/2.E0 
XK3-1.E(K-XK/3.E0 
TOTKE^XK*PRM/6 . E0/G*VP*VP 
GO  TO  <110, 120, 130), IND 

MEAN  PRESSURE  IS  GIVEN 

110  PBRCH-PC*XK2/XK3 
PB-PC/XK3 
RETURN 

BREECH  PRESSURE  IS  GIVEN 

120  PC-PBRCH*XK3/XK2 
PB-PBRCH/XK2 
RETURN 

BASE  PRESSURE  IS  GIVEN 

130  PC-PB*XK3 
PBRCH-PB*XK2 
RETURN 

COMPUTE  SLIP  COEFFICIENT 
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c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


140  BET-1. EO 

IF(ALF.GT. .999E0)  RETURN 

BET-(2 . E0*PRM/C*(PBRCH/PB-1 . EO) -ALF- (PLM+XM0)/C)/(1 . EO-ALF) 
IF(BET.LT.O.EO)  BET-O.EO 
IF(BET.GT.l.EO)  BET-1. EO 
RETURN 

CHAMBRAGE  PRESSURE  GRADIENT  ANALYSIS 

200  IF(KALC.EQ.l)  GO  TO  205 
XP2-XP*XP 
V0L-V0L0+AB*XP 

AJl-AJlO+VOL0/AB*XP+O . 5E0*XP2 
AJ2-V0L*V0L/A?</AB 

AJ3-AJ30+AB*aJ10*XP+VOL0/2 . E0*XP24AB/6 . E0*XP*XP2 
AJ4-AJ4(M-  ( V0L**3 - V0L0**3 ) /3 . EO/AB/AB 
A1T-XK*PRM*AB/G/V0LA0L*  (  AB*VP*VP/VOL*G*AB/PRM*RES  ) 

A2T— XK/AJ2 

BT— XK*PRM/2 . E0/G*VP*VP/AJ2/V0L 

IF(MODE.EQ.l)  GO  TO  240 

XK1-A1T*AJ 1+BT*AJ  2 - A1T*AJ  3/VOL- BT*AJ4/V0L 

XK2-1 . EO - A2T*AJ 1+A2T*AJ  3/VOL 

XK3-1.E0-A2T*AJ1 

XK4— A1T*AJ  1  -  BT*AJ  2 

T0TKE-XK*PRM/2 . E0/G*AJ4/AJ2/VOL*VP*VP 
205  GO  TO  (210, 220, 230). IND 

MEAN  PRESSURE  IS  GIVEN 

210  PBRCH-XK3/XK2*PC+XK3*XK1/XK2+XK4 
PB-(PC+XK1)/XK2 
RETURN 

BREECH  PRESSURE  IS  GIVEN 

220  PC-XK2/XK3*(PBRCH-XK3*XK1/XK2-XK4) 

PB-(PBRCH-XK4)/XK3 

RETURN 

BASE  PRESSURE  IS  GIVEN 

230  PC-XK2*PB-XK1 
PBRCH-XK3*PB-^XK4 
RETURN 

COMPUTE  SLIP  COEFFICIENT 

240  AIT-AIT/XK 
A2T-A2T/XK 
BT-BT/XK 

XK- - ( PBRCH - PB ) / ( A2T*AJ 1*PB+A1T*AJ l+BT*AJ 2 ) 
BET-(PRM/C*XK-ALF- (PLM+XM0)/C)/(1 . EO-ALF) 

IF(BET.LT.O.EO)  BET-O.EO 
IF(BET.GT.l.EO)  BET-1. EO 
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RETURN 

300  WRITE(6,310)  IND 

310  F0RMAT('  GRAD  CALLED  WITH  INVALID  ARGUMENT  IND  -'.IS,' 
*.') 

STOP 

320  WRITE(6,330)  NSTA 

330  F0RMAT(*  GRAD  CALLED  WITH  INVALID  ARGUMENT  NSTA  -'.'IS. 
*'  TERMINATING.') 

STOP 

END 


TERMINATING 
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SUBROUTIHE  SMOOTH 

ROOTINE  TO  FILTER  HIGH  FREQUENCY  COMPCWENTS  OF  DATA  ARRAYS 
COMPLEX  ARRAY. XND 

COMMON/C9/  BALT(1000),BALP(1000).BALPB(1000).MBAL 
DIMENSION  ARG(2000).ARRAY(1000.2).XND(3} 

EQUIVALENCE  (ARG(l) .BALP(l)) 

DATA  XND/(0 . 49965 , 0 . ) . ( -0 . 22227 , 0 . 64240) , ( -0 . 22227 , -0 . 64240)/ 

DO  100  101.2 
IBIAS-1000*(IC-1) 

DO  10  I-l.NBAL 

Il-I+IBIAS 

ARRAY(I.1)-ARG(I1) 

10  CONTINUE 
Jl-2 

J2-IIBAL<1 
DO  30  1-1.3 
K-M0D(I.2)+1 
L-M0D(I+1.2)+1 
DO  20  J-^1.J2 

ARRAY(J  .K)-ARRAY(J .  L)+0 . 5*XND(I)*(ARRAY(J- 1 .  L)-l-ARRAY(J-f  1 .  L)  - 
+  2.*ARRAY(J.L)) 

20  CONTINUE 

ARRAY( 1 . K)-ARRAY( 1 . L) 

ARRAY(NBAL.  K)-ARRAY(NBAL.  L) 

30  CONTINUE 

DO  40  J-l.NBAL 
Jl-J+IBIAS 
BIT<4EAL(ARRAY(  J .  K)  ) 

ARG(J1)-BIT 
40  CONTINUE 
100  CONTINUE 
RETURN 
END 


103 


o  o  o 


SUBROUTINE  REGR 

COMPUTES  LEAST  SQUARES  FIT  TO  (DAT. DATA). 

PARAMETER  (NYI-4.NYPD-4) 

C 

COMMON/C6/  TIME.YI(NYI),YDP(NYPD),YDT(NYI) 

COMMON/C12/  DAT(20) ,DATA(20) .DATEX.AEX.BEX.MDTR.NSTEP 
C 

XNDTR<4n>TR 

TAV-O.EO 

AEX-O.EO 

STA-O.EO 

ST-O.EO 

DO  10  I-1,NDTR 
TAV-TAV+DAT(I) 

AEX-^AEX-fDATA(I) 

STA-STA-»-OAT(  I  )*DATA(  I ) 

ST-ST+DAT ( I ) *DAT ( I ) 

10  CONTINUE 

TAV-TAV/XNDTR 

AEX-AEX/XNDTR 

BEX-( STA-XNDTR*TAV*AEX) / ( ST -XNDTR*TAV*TAV) 
DATEX-AEX-t-BEX*  (  TIME  •  TAV  ) 

RETURN 

END 
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APPENDIX  B 

PHAP:  Code  Listing  and  Description  of  Input 
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Intentionally  left  blank. 
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C*MFST  $ST0RAGE:2 

C  PROGRAM  TO  STUDY  DEPENDENCE  OF  PRESSURE  ON  FRACTION  OF  WORKING 
C  FLUID  CONVERTED  TO  GAS.  ASSUMES  BLAKE  CODE  CHARACTERIZATION  OF 
C  EQUATION  OF  STATE. 

C 

C  VERSION  1.0  (DECEMBER  1989)  WRITTEN  FOR  F77  COMPILER  ON  CRAY. 

C 

C  ***************  *i>0*******ili>i>*»*iliH>A****<!******i>r*it»*«  ***««*  ******** 

C  CODE  MAY  BE  CONVERTED  TO  VERSION  COMPATIBLE  WITH  MICROSOFT  F77 

C  COMPILER  BY  MEANS  OF  GLOBAL  TEXT  EDITING  COMMAND  OF  THE  FORM: 

C 

C  REPLACE  'C*MSFT  '  WITH  *  *. 

C 

C  *****************************"************************************* 

C  MICROSOFT  F77  VERSION  PROMPTS  AT  TERMINAL  FOR  INPUT  DATA  FILE  NAME 

C  (FILNAM)  AND  OUTPUT  FILE  NAME  (OUTFIL) .  THE  OUTPUT  WILL  BE 

C  DIRECTED  TO  THE  PRINTER  IF  OUTFIL  IS  ENTERED  AS  PRN. 

C 

C  '************'***"*************'*'*■*****■**********'********************* 

C  DESCRIPTION  OF  INPUT  DATA  (ASSUMED  TO  BE  LOCATED  IN  FILE  LABELED 
C  FILNAM) : 

C 

C  FILE  1.0: 

C  NVDAT 

C 

C  NCOV 

C 

C  NEOV 

C 
C 

C  FILE  2.0: 

C 
C 

C  VEIN(l) 

C 

C  VEFF(l) 

C 

C  VGAM(l) 

C  VB(1) 

C  VEIN(2) 

C 
C 
C 

C  VB(NVDAT) 

C 

C  NOTES : 

C 
C 
C 
C 


"FILE  COUNTERS"  (315)  ONE  CARD. 

•  NUMBER  OF  SETS  OF  DATA  IN  FILE  2.0. 

NVDAT  MAY  NOT  BE  GREATER  THAN  100. 

•  NUMBER  OF  DATA  IN  FILE  3.0. 

NCOV  MAY  MOT  BE  GREATER  THAN  10. 

•  NUMBER  OF  DATA  IN  FILE  4.0. 

NEOV  MAY  NOT  BE  GREATER  THAN  5. 

"PROPERTIES  OF  FINAL  DECOMPOSITION  OF  WORKING  FLUID 
AS  A  FUNCTION  OF  IXTrAL  INJECTED  PLASMA  ENERGY" 
(4F10.0)  NVDAT  CARDS. 

-  FIRST  VALUE  OF  PLASMA  ENERGY  PER  UNIT  MASS 
OF  WORKING  FLUID  (J/G) . 

-  CORRESPONDING  VALUE  OF  EFFECTIVE  INTERNAL  ENERGY 
(J/G) . 

-  CORRESPONDING  VALUE  OF  RATIO  OF  SPECIFIC  HEATS  ( - ) • 

-  CORRESPONDING  VALUE  OF  COVOLUME  (CC/G) . 

-  SECOND  VALUE  OF  PLASMA  ENERGY  PER  UNIT  NASS 
OF  WORKING  FLUID  (J/G).  STARTS  A  NEW  CARD. 


IF  THE  ENERGY  DENSITY  OF  THE  MIXTURE  OF  PLASMA 
AND  DECOMPOSED  WORKING  FLUID  LIES  OUTSIDE  THE 
RANGE  OF  THE  TABLE  THE  FIRST  OR  LAST  VALUES  ARE 
USED  TO  DETERMINE  THE  EFFECTIVE  PROPERTIES  OF 
THE  MIXTURE. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


THE  ONLY  EXCEPTION  TO  THIS  RULE  OCCURS  FOR  THE 
EFFECTIVE  INTERNAL  ENERGY  AT  VALUES  OF  ENERGY 
DENSITY  LARGER  THAN  VEIN(NVDAT)  .  THE  EFFECTIVE 
INTERNAL  ENERGY  IS  DETERMINED  IN  THIS  CASE  BY 
LINEAR  EXTRAPOLATION  OF  THE  LAST  TWO  VALUES  IN 
THE  TABLE. 

INSIDE  THE  TABLE  RANGE  LINEAR  INTERPOLATI(M(  IS 
USED  TO  DETERMINE  VALUES  OF  ALL  PROPERTIES. 

FILE  3.0:  "VALUES  OF  C/V"  (8P10.0)  ONE  OR  TWO  CARDS. 

DEN(l)  -  FIRST  VALUE  OF  C/V  (G/CC)  .  TH5  RATIO  OF  THE  TOTAL 
MASS  OF  WORKING  FLUID  TO  THE  TOTAL  CHAMBER  VOLUME. 
DEN(2)  -  SECOND  VALUE  OF  C/V. 


DEN(NCOV) 

FILE  4.0:  "VALUES  OF  E/C"  (8F10.0)  ONE  CARD. 

EDEN(l)  -  FIRST  VALUE  OF  E/C  (J/G) ,  THE  RATIO  OF  THE  TOTAL 
ENERGY  TO  THE  TOTAL  HASS  OF  THE  WORKING  FLUID. 
EDEN(2)  -  SECOND  VALUE  OF  E/C  (J/G). 


EOEN(NEOV) 


FILE  5.0:  "PROPERTIES  OF  UNDECC»1P0SED  WORKING  FLUID" 

(3F10.0)  ONE  CARD. 

RHOPO  -  DENSITY  OF  WORKING  FLUID  AT  ONE  ATMOSPHERE  (G/CC) . 

AKl  -  BULK  MODULUS  AT  ONE  ATMOSPHERE  (MPA) . 

AK2  •  DERIVATIVE  OF  MODULUS  WITH  RESPECT  TO  PRESSURE  (•). 

***************************************************************** 
CHARACTER* 10  FILNAM , OUTFIL 

COMMON/Cl/  RHOPO, AKl, AK2 


DIMENSION  PRES(50, 10) ,DEN(10) ,EDEN(5) ,VEIN(100) ,VEFF(100) , 
*  VGAM(100),VB(100) 

C 

C  NAME  FILES  FOR  INPUT  AND  OUTPUT  ON  MICROCOMPUTER 
C 


C*MFST 

WRITE(*,1) 

C*MFST 

1 

FORMAT ('  ENTER  INPUT  DATA  FILE  NAME: 

C*MFST 

READ(*,2)  FILNAM 

C*MFST 

2 

FORMAT(AIO) 

C*MFST 

OPEN( 5 , FILE-FIINAM) 

C*MFST 

WRITE(*,3) 

C*MFST 

3 

FORMAT('  ENTER  OUTPUT  FILE  NAME:  '\) 

C*MFST 

READ(*,2)  OUTFIL 

C*MFST 

OPEN(6 , FILE-OUTFIL) 

C 
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C  READ  AND  WRITE  INPUT  DATA 
C 

READ(5,100l)  NVDAT.NCOV.NEOV 

READ(5, 1005)(VEIN(I) ,VEFF(I) ,VGAM(I) .VB(I) , I-l.NVDAT) 

READ(5,1002)  (DEN(I) , I-l ,NCOV) 

READ(5.1002)  (EDEN(I) . I-l .NEOV) 

READ(5.1002)  RHOPO.AK1.AK2 
WRITE(6,2120)  NVDAT.NCOV.NEOV 

WRITE(6.2130)  (VEIN(I) .VEFF(I) .VGAM(I) .VB(I) , I-l.NVDAT) 

WRITE(6 . 2132)  (DEN(I) . I-l .NCOV) 

WRITE(6.2134)  (EDEN(I) .1-1. NEOV) 

WRITE(6.2060)  RHOP0.AKl.AK2 
C 

C  CALCULATE  VALUES  OF  PRESSURE 
C 

DO  300  K-l.NEOV 
EOC-EDEN(K) 

DO  200  J-l.NCOV 
COV-DEN(J) 

DO  100  1-1.50 
ALF-0.02E0*  FLOAT(I) 

PLEIN-EOC/ALF 

CALL  INTERP ( PLEIN , EMIX . DXEMIX . VEIN . VEFF , 1 . NVDAT) 

CALL  INTERP ( PLEIN . GAM . DXGAM . VEIN . VGAM . 0 . NVDAT ) 

CALL  INTERP ( PLEIN. B.DXB. VEIN. VB.O. NVDAT) 

BIT1-(GAM- 1 . EO)*ALF*EMIX 
BIT2-1 . EO/COV-B*ALF 
BIT3— (l.EO-ALF) 

PC-O.EO 

BITA-BIT2+BIT3/RHOP0 
IF(BIT4.LT.l.E-3)  GO  TO  90 
PC-BIT1/BIT4 
NIT-O 

30  CALL  EOSLIQ(PC.RHOP.DRHOP.EP.DEP) 

BIT4-BIT2+BIT3/RHOP 

C 

C  PRESSURE  AND  LIQUID  DENSITY  ARE  COUPLED  NONLINEARLY 
C 

50  FP-PC-BIT1/BIT4 

IF(  ABS(FP).LT.1.E-4*PC)  GO  TO  90 
IF(NIT.LT.50)  GO  TO  80 
WRITE(6.60;  FP.PC 

60  FORMAT('  MORE  THAN  50  ITERATIONS  TO  DETERMINE  PRESSURE  IN  DEPEND. 
★TERMINATING, '/ 

*'  FP  -',G15.6,'  PC  -’,G15.6) 

STOP 

80  DFP-1 . EO-DRHOP*BIT1*BIT3/(BIT2*RHOP+BIT3)**2 
PC-PC-FP/DFP 
NIT-NIT+1 
GO  TO  30 
90  PRES(I,J)-PC 
100  CONTINUE 
200  CONTINUE 
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PRINT  TABLE  OF  RESULTS 
WRITE(6,2090) 

WRITE(6,3000)  EDEN(K) , (DEN(J) .J-l.NCOV) 

DO  260  1-1,50 
ALF-0.02EO*  FLOAT(I) 

WRITE(6 , 3010)  ALF, (PRES(I , J) , J-l,NCOV) 

260  CONTINUE 
300  CONTINUE 
STOP 

FORMAT  STATEMENTS 

1001  FORMAT(16I5) 

1002  FORMAT (8F10.0) 

1005  FORMAT(4F10.0) 

2060  F0RMAT('  '// 

*'  ’ .5X, 'PROPERTIES  OF  UNDECOMPOSED  WORKING  FLUID'// 

*'  DENSITY  OF  FLUID  AT  1  ATM.(G/CC)  '  FIO  4/ 

*'  BULK  MODULUS  AT  1  ATM. (MPA)  '[fIO.1/ 

*'  DERIVATIVE  OF  MODULUS  W.R.T.  PRESSURE(-)  ',’f10.3//) 

2120  PORMAT('  PRESSURE  AS  A  FUNCTION  OF  lUACTION  OF  WORKING  FLUID  CONVE 
*RTED  TO  GAS'// 

*'  NUMBER  OF  DATA  TO  DETERMINE  THERMOCHEMISTRY  ',15/ 

*'  NUMBER  OF  VALUES  OF  C/V  ',15/ 

*'  NUMBER  OF  VALUES  OF  E/C  M5//) 

2130  PORMATC'  '  ,5X, 'ENERGY(J/G)  '  ,5X, 'EFF.ENGY.  (J/G) '  ,8X, 'GAMMA(-)  ' , 

*5X. ' COVOLUME(CC/G) ' //< '  ' . 2F15 . I , 5X, 2F15 . 4) ) 

2132  FORMATC  V/‘  VALUES  OF  C/V  (G/CC) :'//('  ',10G12.5)) 

2134  F0RMAT('  '//'  VALUES  OF  E/C  (J/G) :'//('  ',5G12.5)) 

2090  FORMAT('l') 

3000  FORMAT('  ',  lOX ,' PRESSURE  AS  A  FUNCTION  OF  ALPHA  AND  C/V  FOR  E/C 
*,G12.5//'  ' ,6X, 'ALPHA' ,40X, 'C/V'/'  ' ,10X,10F10.3//) 

3010  P0RMAT('  ' ,F10.4,10F10.1) 

END 
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SUBROUTINE  1NTERP(X,  Y,DYDX.XA,YA.NXTRAP,N) 

C 

C  INTERPOLATE  (XA.YA)  ARRAYS  TO  GET  Y(X)  AND  DERIVATIVE  DYDX. 

C  IF  X  IS  OUTSIDE  RANGE  OF  XA.  Y  IS  SET  EQUAL  TO  FIRST  OR  LAST 
C  VALUE  OF  YA  UNLESS  NXTRAP  -  1  IN  WHICH  CASE  LINEAR  EXTRAPOLATION 
C  IS  USED  FOR  X  >  XA(N) . 

C  IF  N  IS  LESS  THAN  TWO,  Y  IS  SET  EQUAL  TO  ZERO. 

C 

DIMENSION  XA(N) .YA(N) 

C 

Y-O.EO 
DYDX-O.EO 

IF(N.LT.2)  GO  TO  50 
IF(X.LT.XA(1))  GO  TO  30 
IF(X.GT.XA(N))  GO  TO  40 
DO  10  1-2. N 

IF(XA(I).GT.X)  GO  TO  20 
10  CONTINUE 
GO  TO  40 

X  IS  IN  TABLE  RANGE.  INTERPOLATE. 

20  BIT-XA(I)-XA(I-1) 

IF(BIT.GT.1,E-10)  DYDX-(YA(I)-YA(I-1))/BIT 
Y-YA(I - 1)+(X-XA( I - 1 ) )*DYDX 
GO  TO  50 

X  IS  OUTSIDE  TABLE  RANGE.  USE  FIRST  OR  LAST  VALUE  OF  YA. 

30  Y-YA(l) 

GO  TO  50 

40  IF(NXTRAP.EQ.l)  GO  TO  45 
Y-YA(N) 

GO  TO  50 
45  I-N 

GO  TO  20 
50  RETURN 
END 
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SUBROUTINE  EOSLIQ( PIN . RHOP . DRHOP . EP . DEP) 


DETERMINES  DENSITY  (RHOP)  AND  INTERNAL  ENERGY  (EP)  OF  WORKING 
FLUID  AS  A  FUNCTION  OF  PRESSURE  (PIN)  TOGETHER  WITH  DERIVATIVES 
DRHOP  AND  DEP.  PMAP  VERSION  SETS  EP-DEP-O.EO 

COMNON/Cl/  RHOPO.AK1.AK2 

P-PIN-O.IOIEO 

IF(P.GT.O.EO)  GO  TO  10 

RHOP-RHOPO 

DRHOP-O.EO 

EP-O.EO 

DEP-K).E0 

RETURN 

10  IF(  ABS(AK2).GT.1.E*6)  GO  TO  20 
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Nominal  Data  Base  for  Unped  Parameter  Study 
of  Ballistic  Implications  of  Rate  of  Mixing 
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TEST  CASE  BASED  ON  40HTC3 

TIHE(HS)  P-BR(HPA)  P-BASE(HPA)  TEHP(K)  YFRC-)  PR.  DISP(CM)  PR.  VEL(H/S)  ALPHA(-)  BETA(-)  E-mC(-) 


oooooooooooooooooooooooooooooooooooooooooooooeoooo 

oooooooooooooooooooooooooooooooooooooooooooooooooo 

OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO 


oooooooooooooooooooooooooooooooooooooooooooooooooo 

oooooooooooooooooooooooooooooooooooooooooooooooooo 

oooooooooooooooooooooooooooooooooooooooooooooooooo 


aoflo«tt«a<o«oas«««0«o«o«0«B«a«o««oaaM«»«o«»flo«oao«flQ«»M«»aoo^^0^0A^94a^o^«B^<^m»«B0Ba«^^ 


r*  4^  ^  ^ 

«»  49  tfl  ^  ^ 

c«  c#  Cl  r«  n 


o 

o 

o 


o 

o 

o 


o  o 
o  o 
o  o 


o  o  o 
o  o  o 
o  o  o 


o  o 
o  o 
o  o 


o  o  o 
o  o  o 
o  o  o 


o  o  o 
o  o  o 
o  o  o 


oooooooooo 

oooooooooo 

oooooooooo 


o  o  o 
o  o  o 
o  o  o 


ooooooooooooooooeoooo 

ooooooooooooooooooooo 

ooooooooooooooooooooo 


O  9»  40  flO 

r»  ^  C4 


cs  n  r<«  r« 


AOAO«o^«nci4>49t 

O0^40r*«oin^c4 

^ooooooo 

cicic«r«r«c4r4c^ 


A0^«n«-iO4or»«n^ 
9«49r>^w>«nr4^0 
^  ^  9%  ^  Ob  ^ 


tfB  Ok  <*B  r*  o 

O  Ok  m  40  « 


r*  o  tn  BO  ^  *>4 


iB0O^r-<^B0O404 

)40<o«nO40tn«nof 


>C4r4r4r4r«C4r<4c«c«r4« 


ir*r-r>r»r*r»r«BOBOBO«OBOBO 


r*4^^40tOi/>4iPmc«*i4O0k0B0k0kOk^OO«^r4«^0>ro0k^p*B^^«ncBr4Or>^v4r^r440Ci^O9r<>^r4^«nr>40 

^4O^lAOBr>r>«>4lO9B(n«OO^4OC«BO«^«f^OB«^r>*«H«B>OB^0r4^rO^^fnOlBO*nOk<OC4OB«nC44O^Or*>4^^«A»^ 


I  ^  ^  o  o 


wi«r)^^<ncir«^^o^^40 

r*r*r»r*r*t-r»c»r«‘r*BO'OBO 


«n  in  lA  4P 
lO  ^  ^  >o 


«nmC4€^C4^^000 


oooooooo 

Oo4r«4n9«n«or> 

tn«nintnin*ninin 


o  o 

40  OB 

«n  in 


o  o  o  o  o  o 

O  c«  «»B  ^  in 
lO  *0  «0  BO  40  «o 


o  o  o  o  o  e 
^  r**  40  OB  o 
BO  BO  bo  bo  r* 


oooooooooooooooooooooooooo 

9O»Bor>‘40^O*^c4««B^in«r»400BO«-4Ci«*BB0inB0r»400B 

r^r*r"r*>P*r>4O4O4O4O4O4O4O4O4D4O^OB0B^OBOB^^^OB 


<*4<^«r4C4r4C4C«r4<Mc«C4r4c«r4ncic«r4C«dr>«r4r4Cic«r4r4dc4r4r«ciC4r4r4i^r4r4r«cic*4C4r4r4C4r404C«r4r4 


129 


TEST  CMI  BASKO  ON  SONTCS 
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APPENDIX  E 

Nonlnal  DaCa  Base  for  XXTC  Sinulatlon  of  ETC 
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